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ABSTRACT 


The  subject  of  this  investigation  is  the  deformation  and  consequent  lack  of  aim 
that  a  gun  tube  undergoes  as  a  result  of  its  initial  imperfections.  By  initial  imperfections, 
we  mean  the  non-straightness  of  the  centroidal  axis  of  the  gun  tube.  In  reality,  the 
centroidal  axis  is  a  general  space  curve.  For  the  sake  of  analysis,  however,  it  is  assumed  to 
be  helical  in  shape.  The  gun  tube  is  modeled  as  a  rod  of  linearly  elastic  material,  encased  in 
a  perfectly  straight  cylindrical  thermal  shroud  to  protect  it  from  non-uniform  temperature 
changes.  The  interaction  between  these  two  elements,  namely,  the  shroud  and  the  gun 
tube  renders  the  problem  of  finding  the  muzzle  end  displacement  of  the  gun  tube  statically 
indeterminate. 


Two  soLcions  to  the  problem  are  proposed  herein,  using  different  approaches.  In 
the  first,  a  more  exact  solution  for  the  gun  tube  end  displacement  is  obtained  by 
formulating  an  expression  for  the  complementary  strain  energy  of  deformation  of  the  rod 
and  then  employing  Castigliano's  Second  Theorem,  or  as  it  is  better  known,  by  the  method 
of  minimum  complementary  strain  energy.  In  the  second  approach,  the  differential 
equations  for  the  problem  formulated  by  Kingsbury  are  employed  and  a  general  finite 
element  of  variable  curvature  and  torsion  is  developed.  This  element  is  then  used  to 
compute  a  set  of  displacements  for  independent  confirmation.  The  element  formulated 
may  be  used  to  analyze  tubes  or  rods  of  arbitrarily  varying  cross-section,  curvature  and 
torsion.  The  report  includes  surface  plots  from  a  pwametric  study  of  the  end  displacement 
as  a  function  of  initial  curvature  and  torsion  of  the  gun  tube  and  the  conclusions  drawn 
therefrom.  Code  for  a  finite  element  program  that  uses  the  element  developed  has  been 
provided  for  implementation  by  the  reader. 


CHAPTER  1 


INTRODUCTION 


1.1  Defining  the  Problem 


The  work  described  in  the  pages  to  follow  is  aimed  at  modeling  the  imperfections 
in  the  straightness  of  gun  tubes.  Gun  tubes  are  often  subjected  to  large  variations  in 
temperature  which  arise  from  several  sources,  for  example,  heating  due  to  the  firing  of 
rounds  of  ammunition,  cooling  due  to  rainfall  or  wind,  and  the  effect  of  field  temperature. 
For  this  reason,  they  are  encased  in  a  protective  outer  cylinder  called  a  thermal  shroud, 
whose  function  is  to  distribute  temperature  changes  uniformly. 

The  gun  tube  and  its  thermal  shroud  are  attached  to  a  common  structural  base 
near  the  receiver  end  of  the  gun  and  are  again  structurally  connected  near  the  muzzle.  The 
latter  connection  allows  relative  motion  of  the  end  of  the  shroud  to  be  transmitted  to  and  to 
be  constrained  by  the  enclosed  gun  tube.  Since  the  temperature  changes  are  distributed 
uniformly  around  any  circumference,  the  shroud  does  not  undergo  bending  deformations 
due  to  thermad  strains  but  expands  or  contracts  axially. 

When  the  gun  tube  and  shroud  assembly  is  subject  to  a  change  in  temperature, 
say  an  increase,  the  shroud  undergoes  thermal  expansion  and  thus  exerts  a  force  on  the  gun 
tube.  If  the  gun  tube  were  to  be  perfectly  strught,  this  would  merely  cause  it  to  extend 
along  its  centroidal  axis,  with  no  change  in  its  um.  However,  in  reality,  no  gim  tube  is 
perfectly  strught  so  that  a  change  in  temperature  of  the  shroud  may  produce  not  only  an 
extension  of  the  tube  in  the  direction  of  its  length  but  also  bending  deformation,  which 
affects  its  aim.  We  endeavor,  then,  to  find  the  displacement  and  rotation  of  the  end  of  the 
gun  tube  as  a  function  of  the  imposed  temperature  change. 
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2 


1.2  The  Analytical  Model 

The  centroideil  axis  of  the  gun  tube  would,  in  reality,  be  an  arbitrary  space  curve. 
For  the  purpose  of  analysis,  it  is  assumed  to  be  helical  in  shape.  It  is  also  assumed  that  the 
cross-section  of  the  helical  rod  is  perfectly  circular,  so  that  the  effects  of  out-of-roundness  or 
twist  of  the  tube  are  not  considered. 

The  gun  tube,  then,  is  modeled  as  a  helical  rod  of  linearly  elastic  material,  with  a 
constant  circular  cross-section,  which  is  contained  within  a  thin-walled  right  circular 
cylinder  representing  the  thermal  shroud.  Both  the  rod  and  the  cylinder  are  considered 
fixed  at  one  end  to  a  massive  common  btise.  The  conunon  attachment  at  the  other  end  of 
the  rod  is  such  that  the  relative  motion  of  the  end  of  the  cylinder  due  to  thermal  strain,  can 
be  imputed  to  -  and  is  constrained  by  -  the  rod  at  their  point  of  intersection.  Two  cases  of 
compatibility  of  gun  tube  and  shroud  end  displacements  are  considered,  the  first  in  which 
the  shroud  is  assumed  to  have  zero  lateral  stiffness,  and  the  second  in  which  the  shroud  has 
a  finite  non-zero  lateral  stiffness.  In  the  former  case,  the  gun  tube  end  displacements  and 
rotations  are  not  affected  by  the  transverse  rigidity  of  the  cylinder.  In  both  cases,  however, 
these  displacements  and  rotations  are  independent  of  the  torsional  rigidity  of  the  cylinder. 

1.8  Methods  of  Solution 

Two  solutions  to  the  problem  are  proposed  herein,  using  different  approaches.  In 
the  first,  more  exact,  solution,  the  end  displacement  is  obtained  by  formulating  an 
expression  for  the  strain  energy  of  deformation  of  the  rod  in  terms  of  the  components  of  the 
applied  force  and  moment  and  then  employing  Castigliano’s  Second  Theorem,  or  as  it  is 
better  known,  by  the  method  of  minimum  complementary  stridn  energy.  In  the  second 
approach,  the  differential  equations  for  the  problem  formulated  by  Kingsbury  [Kingsbury 
84]  are  employed  and  a  general  finite  element  of  variable  curvature  and  ttwsion  is 
formulated  using  the  variational  method.  This  element  is  then  used  to  independently 
confirm  the  results  of  the  analytical  model. 
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1.3.1  An  Exact  Solution 

The  required  displacements  are  obtuned  most  ea^y  by  employing  an  energy 
method;  to  be  specific,  by  using  Castigliano’s  Second  Theorem  which  is  stated  as  follows:  If 
a  linearly  elastic  structure  is  subjected  to  a  set  of  loads,  the  displacement  of  any  load  in  its 
direction  is  equal  to  the  partial  derivative  of  the  complementary  strain  energy  with  respect 
to  that  load.  In  order  to  use  Castigliano’s  theorem,  an  expression  for  the  complementary 
str2un  energy  of  the  rod  must  be  formulated.  This  is  done  as  follows:  First,  the  helix 
representing  the  centroidal  axis  of  the  gun  tube  is  parameterized  in  terms  of  the  arc  length 
and  expressions  are  obtained  for  the  vector  tangent,  binormal  emd  principal  normal  to  the 
curve  at  any  point  along  its  length.  Second,  expressions  for  the  intem^ll  force  and  moment 
components  along  the  local  tangent,  binormal  wd  principal  normal  that  result  from 
applying  a  force  and  a  moment  to  the  end  of  the  rod  are  obtained  in  terms  of  the  applied 
force  and  moment.  Thirdly,  the  complementary  strain  energy  of  the  rod  is  written  in  terms 
of  the  internal  force  and  moment  components,  and  therefore,  in  terms  of  the  applied  force 
and  moment.  The  required  scalar  displacements  are  finally  obtained  by  differentiating  the 
complementary  strain  energy  with  respect  to  each  component  of  the  applied  force  and 
moment. 


Since  the  problem  is  statically  indeterminate,  the  applied  force  is  not  known  but  is 
found  by  enforcing  the  compatibility  of  gun  tube  and  shroud  end  displacements.  This  force 
is  then  used  to  compute  the  required  displacements.  A  parametric  study  of  the  force  and 
displacement  components  is  performed  in  which  each  component  is  plotted  as  a  function  of 
the  initial  curvature  and  torsion  of  the  helical  rod.  This  is  done  for  two  cases  of 
compatibility  of  displtu:ements  of  the  gun  tube  and  the  thermal  shroud  at  the  muzzle  end; 
one  in  which  the  thermal  shroud  is  assumed  to  have  only  axial  stifhess  and  no  transverse  or 
rotational  stifhess  and  the  other  in  which  the  shroud  has  a  finite  non>sero  transverse 
sti&ess.  Chapters  2  and  3  present  the  formulation  of  the  striun  energy  and  equilibrium 
equations  while  Chapter  4  deals  with  the  imposition  of  compatibility  of  displacements  of 
the  muzzle  ends  of  the  gun  tube  and  the  thermal  shroud,  and  consequent  resolution  of 
statical  indeterminacy.  Chapter  5  contains  the  above>mentioned  parametric  study. 
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1.8.2  A  Finite  Element  Solution 

In  Chapter  6,  an  alternative  solution  is  presented  which  involves  the  development 
of  a  finite  beam  element  whicll  is  permitted  to  have  variable  curvature  and  torsion  over  its 
length.  However,  the  rate  of  variation  of  the  curvature  and  torsion  is  assumed  to  be 
constant  in  each  element.  By  the  use  of  an  assemblage  of  such  finite  elements,  rods  of 
arbitrarily  varying  space  curvatures  can  be  analyzed. 


The  element  is  formulated  using  the  variational  method.  The  equations  of  motion 
of  a  space-curved  rod  have  been  formulated  by  Kingsbury  [Kingsbury  84]  as  a  set  of  four 
coupled  differential  equations  in  four  unknown  displacements.  Every  point  on  the  rod  has 
six  degrees  of  freedom:  three  translations  and  three  rotations.  Two  of  these  are  expressed  in 
terms  of  the  four  unknowns  appearing  in  the  differential  equations. 

The  development  proceeds  as  follows:  First,  a  polynomial  displacement  field  is 
assumed,  which  satisfies  the  differential  equations  exactly  in  the  case  in  which  curvature 
and  torsion  are  identically  zero,  that  is,  in  the  case  of  a  straight  rod.  Next,  the  constant 
coefficients  appearing  in  the  assumed  polynomials  are  found  in  terms  of  the  nodal  values  of 
displacements,  and  these  polynomials  are  rewritten  to  yield  the  shape  functions.  Then,  the 
expression  for  linear  strain  energy  density  of  the  rod  developed  by  Tsay  amd 
Kingsbury  [Tsay  86]  is  employed  as  follows.  The  displacements  and  their  derivatives 
appealing  in  this  expression  are  substituted  with  their  equivalents  in  terms  of  the  assumed 
displacement  field  and  the  strain  energy  density  is  integrated  over  the  length  of  the  rod  to 
yield  the  strain  energy  of  the  rod.  Finadly,  the  strain  energy  is  differentiated  with  respect  to 
the  nodail  values  of  the  displacements  to  obtain  the  elements  of  the  stifhess  matrix. 

Lastly,  the  finite  element  is  used  to  independently  confirm  the  values  for 
displacements  obtained  from  the  first  method.  Code  for  a  finite  element  program  has  also 
been  provided  in  Appendix  D  for  implementation  by  the  reader. 


CHAPTER  2 


INTERNAL  FORCES  IN  THE  HELICAL  ROD 


2.1  Parametric  Representation  of  the  Helix:  Vector  Tangent,  Binomoal  and 
Principal  Normal 


2.1.1  Parametric  Representation 

The  geometry  of  the  helical  rod  is  shown  in  Figure  2.1.  Its  shape  is  governed  by 
two  parameters:  a,  which  is  the  length  of  the  circular  arc  representing  the  projection  of  the 
helical  centroidal  axis  of  the  rod  on  the  XY  plane,  and  the  angle  subtended  by  this  arc  at 
the  origin.  The  points  on  the  helix  can  be  described  in  terms  of  the  arc  length,  s,  by  a 
position  vector,  r,  given  by: 


We  also  have  the  following  expressions  for  the  first  and  second  derivatives  of  r 
with  respect  to  s: 


r  =  —  =  - -sm(— )i -cos(— )j  + -k 
da  L  L  l’  L 


r' ' 


(2.2) 


(2.3) 
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(i)  The  "Straightened"  Helical  Arc 
Z 


(ii)  The  Helical  Arc  and  the  Global  Coordinate  System 


Figure  2-1:  Geometry  of  the  Helical  Rod 
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2.1.2  Vector  Tangent,  Binormal  and  Principal  Normal 


These  three  characteristic  vectors  of  a  space  curve  form  what  is  sometimes  referred 
to  as  the  trihedral  associated  with  the  curve  [SokolnikofF  661.  They  are  shown  in  Figure  2.2, 
eind  at  any  point  along  the  curve,  they  are  given  in  terms  of  the  derivatives  of  the  position 
vector,  r,  by  the  following  relations: 


Tangent  r. 


?»)  = 

r  (s) 


Binormal  i3: 


t'(s)x  t'’{s) 
r'(s)x  r"(3); 


Principal  Normal  u: 
i/{s)  =  $ls)  X  r{s) 


(2.4) 


(2.5) 


(2.6) 


Substituting  the  derivatives  in  the  above  expressions,  we  find  the  following 
expressions  for  the  three  unit  vectors; 

r(j)  =  -  -  sin  (— )  i -t- -  cos  (— )  j  +  -  k  (2.7) 

'  L  L  L  L 

0{s)  =  -sin(— )i  -  -cos(^)j  + -k  (2.8) 

L  l'  L  L 

v[»)  =  -  co8(^)i-  sin(^)j  (2.9) 


The  vector  tangent,  binormal  and  principal  normal  form  a  local  coordinate  system 
which  varies  with  every  point  on  the  curve.  We  assume  that  the  tangent  points  in  the  local 


X(S-i^ds) 


Figure  2-2:  Vector  Tangent,  Binormal  and  Principal  Normal 
z  or  ;  direction.  Therefore,  the  principal  normal  and  the  binomoal  point  in  the  local  x  (() 
and  y  {tj)  directions  respectively.  In  the  next  section,  we  shall  use  these  three  unit  vectors 
to  evaluate  the  components  of  the  internal  force  and  moment  along  the  directions  of  the 
local  rf  and  ^  axes. 

2.2  Internal  Force  and  Moment  Components  at  a  Generic  Point 

As  mentioned  in  Chapter  1,  the  muzzle  end  of  the  gun  tube  is  structurally 
connected  to  the  thermal  shroud  in  a  manner  that  allows  motion  of  the  end  of  the  shield  to 
be  transmitted  to  and  also  to  be  construed  by  it.  Also,  the  receiver  end  of  the  gun  and 
the  shroud  are  connected  to  a  common  masnve  base.  In  other  words,  for  our  analyra,  we 
need  to  find  the  displacement  of  the  firee  end  of  a  helical  rod,  the  other  end  being 
cantilevered.  Having  done  this,  we  may  then  apply  compatibility  of  displacements  of  the 
shroud  and  gun  tube  and  find  the  force  that  results  fium  their  interaction.  This 
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corresponds  to  the  applied  force  F  in  the  discussion  to  follow. 


2.2.1  Internal  Force  and  Moment 


We  now  evaluate  the  internd  force  and  moment  components  at  a  generic  point. 
Figure  2.3  shows  the  forces  acting  on  the  &ee  end  of  the  rod  and  the  corresponding 
reactions  at  the  fixed  end  as  represented  in  the  global  XYZ  system. 

The  applied  force  at  point  B,  that  is,  at  s  =  L,  is  F  given  by: 

F  =  (2.10) 

The  applied  moment  at  point  B  is  M  given  by: 

M  =  Mxi-Myj  +  Mzk  (2.11) 

The  reactions  at  the  fixed  end  A,  that  is,  at  s  =  0,  are  a  force  F^,  and  a  moment 
F«  is  found  from  force  equilibrium  as: 

Ffi  =  “  F  =  -Fxi-Fyi-F^k  (2.12) 

is  found  from  moment  equilibrium  by  summing  up  moments  about  the  point  A,  and  is 
given  as: 


'M.ji  =  -  M-  A  X  F,  where  A  =  r(L)  -  r(0), 

=  [-  Mx+Fyb-^^t^]i-  [J\fy+Fxfc  +  ^(l-cos^)]j 

<t>  <t> 

+  [-  (2.13) 

<b  <j> 


Consider  the  free  body  shown  in  Figure  2.4.  Now,  from  force  equilibrium,  we  have 
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z 


Reiction  force 
F«-Fxi-Fyj-F2k 


Figure  2-8:  Forces  on  the  Helical  Rod 


Z 


Resction  force  *<•>  “ 

F--F«|.Fyj.Fxk 


Figure  2-4:  Free  Body  Diagram 
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the  interned  force  at  a  generic  point  P,  ,  in  terms  of  the  eurc  length  along  the  rod,  5,  as: 

-  ~^R  =  + (2-14) 

Setting  the  sum  of  moments  about  the  point  A  to  zero  (moment  equilibrium),  we 
find  the  interned  moment  at  point  P,  ,  to  be  given  as: 

M^(s)  =  -M;j-Ar(s)xF^(s)  (2.15) 

where  Ar(3)  is  the  lever  eirm  given  by: 

Ar(s)  =  r(3)-r(0)  (2.16) 

The  complete  expression  for  the  internal  moment  acting  at  the  generic  point  P,  in  terms  of 
the  arc  length,  s,  is: 

Mrt(«)  =  .  Mx~  ^  ~  «)  -—  ( sin -  sin  (^)  ]  i 

L  <p  L 

+  s)-  ^(cos(i>  -  cos(^))]} 

L  <p  L 

-r  M2  ~  ( sin  ^  -  sin  (— ) )  -h  ( cos  ^  -  cos  (^) )  ]  k  (2-17) 

<t>  L  L 

2.2.2  Local  Components  of  Internal  Force  and  Moment 

The  components  of  the  internal  force  and  moment,  Vj,  V,^,  and  Mj,  M,,, 
respectively,  may  be  found  by  taking  the  dot  product  of  the  vector  in  question  with  the  unit 
vectors  in  the  local  coordinate  system,  namely,  the  vectors  v,  0  and  r.  Therefore,  we 
have  the  following  scalar  equations  for  the  internal  force  and  moment  components  as 
functions  of  the  arc  length  s: 

V((s)  =  -  FxCM  (y )  -  Fysin  (^) 


(2.18) 


12 


Vn(«)  =  Sm  (— )  -  — £-  COS  (— )  +  — ^ 

^  L  ^l'  L  ^l’  L 

V(,)  =  -!x±An<t)+EyIc„ip)+!^ 

^  L  l’  L  l'  L 

M({s)  =  [-  Mx^^{L-  3)--^{sm<^-sin(^)}]co8(^) 
L  0  L  L 


-  -  M. 


-£^{L-  3)+:^{cos.^-co8(^)}]sm(^) 


<t> 


Mr,{s)  =  [Mx~  3)+^{sin<^-sui(^)}lisiii(^) 

L  <p  L  L  L 


-  My-  {L-  3)-!-^^{cOS<^-  C08(^)}  ]-^C08{^) 

L  4>  L  L  L 


•r  '  sm<^-  8in(^)}  -j-^^{  co8<^-  cos(^)}  \  ^ 

<p  L  4>  L  L 

M^{s)  =  ;-  Mx  +  ^{L-  3)-:^{9m(^-  8m(^)}  lysinA 

L  0  L  If  h 


\My  +  ^^{L-  3)-^^{  CO8  0-  COs(^)}  ]^C09  (^) 

L  <t>  L  L  L 


[M^-  ^^{sin<t>-  siii(^)}  co8«^-  co9(^)}  l| 

0  L  4>  L  L 


(2.19) 

(2.20) 


(2.21) 


(2.22) 


(2.23) 


In  the  following  chapter,  we  find  an  expression  for  the  complementary  strain 
energy  of  the  deformed  helical  rod  in  terms  of  the  above  components  of  the  internal  force 
and  moment  acting  at  a  croes-section,  and  then  apply  Castigliano’s  Second  Theorem  to  find 
the  displacements  at  the  free  end  of  the  helical  rod,  and  thereby,  the  elements  of  the 
corresponding  flexibility  matrix. 


CHAPTER  S 


THE  STRAIN  ENERGY  OF  THE  ROD  AND  ITS  END  DISPLACEMENTS 


3.1  Strain  Energy  of  the  Deformed  Helical  Rod 


3.1.1  Components  of  Stress 


The  components  of  the  Internal  force  and  moment  at  any  cross-section  are  defined 
in  terms  of  the  local  stress  components  by  the  following  integrals  evaluated  over  the  area  of 
that  cross-section: 


=  - 

Ja 

/"  (3-1) 

Even  though  we  do  not  have  relations  for  the  components  of  stress,  viz.  and  in 

terms  of  the  components  of  the  internal  force  and  moment  at  a  croes-section,  we  assume 
relations  (analogous  to  known  strength  of  materials  relations)  which  are  such  that  on 
substitution  in  the  defining  formulae,  they  yield  identities.  The  assumed  formulae  for  the 
stress  components  are: 
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VM,  ,  Vf 


3.1.2  Strain  Energy  in  Terms  of  Internal  Force  and  Moment  Components 

We  now  derive  an  expression  for  the  complementary  str^^n  energy  of  the  deformed 
helical  rod  in  terms  of  the  local  components  of  the  interned  force  and  moment  acting  at  a 
cross-section.  Assuming  the  material  constituting  the  rod  to  be  linearly  elastic,  the  strain 
energy  of  the  rod  is: 


=  -  / 
2  Iv 


dV 


The  stress  and  strain  vectors  in  the  above  expression  are: 

<r  ^  =  :  j 

«  ^  =  >{f  S'/  Sf  Hn  Sf  Sf  1  (3.4) 

As  per  the  assumptions  made  in  the  analysis  of  a  general  curved  and  twisted 
rod  Kingsbury  84 we  have: 


-  S'?  -  Un  -  ^ 


With  these  assumptions,  the  applicable  stress-strain  relations  are: 


Sf 


(3.6) 
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Substituting  from  equations  (3.4),  (3.5)  and  (3.6)  into  Equation  (3.3),  we  get  the  following 
expression  for  the  complementary  strain  energy; 


Since  the  bar  is  of  constant  circular  cross-section  of  area  A,  we  have  the  following 
section  properties: 


hr,  =  0 


j  ndA  =  0 

A  A 


(3.8) 


Substituting  equations  (3.2)  and  (3.8)  into  Equation  (3.7),  we  arrive  at  the 
following  expression  for  the  complementary  strain  energy  of  the  deformed  helix  m  terms  of 
the  components  of  the  internal  force  and  moment  represented  in  the  local  coordinate 
system: 


(3.9) 


It  may  be  noted  that  the  contributions  to  the  shear  stress  components  in 
equations  (3.2)  from  the  local  force  components  and  are  assumed  to  be  constant  over 
the  cross-section.  More  accurate  results  may  be  obtained  by  assuming  a  parabolic 
distribution  of  shear  stress  instead  of  an  average  value. 
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S.l.S  Strain  Energy  in  TemoB  of  Applied  Force  and  Moment  Components 

In  Chapter  2,  we  found  expressions  relating  the  local  components  of  the 
internal  force  and  moment  at  a  cross-section  to  the  global  (XYZ)  components  of  the  applied 
force  F  and  moment  M.  We  now  combine  those  relations,  namely,  equations  (2.18)  - 
(2.23),  with  Equation  (3.9),  which  is  the  expression  for  the  complementary  strain  energy  of 
the  rod  in  terms  of  the  components  of  the  internal  force  and  moment  at  a  cross-section. 
Integrating  the  expression  so  obtained  over  the  length  of  the  rod,  L,  we  find  the 
complement2u:y  strain  energy  of  the  deformed  rod  in  terms  of  the  geometrical  parameters, 
viz.  a  and  0,  and  the  XYZ  components  of  the  applied  force  and  moment.  (Note  that  b  is 

not  an  independent  parameter,  but  is  defined  in  terms  of  a,  as  \J  1?  -  a,  L  being  held 

constant.)  The  complete  expression  for  the  complementary  str^  energy  of  the  rod  may  be 
found  in  Appendix  A.  We  are  now  equipped  to  ctdculate  the  displacements  that  occur  at 
the  free  end  of  the  rod. 

3.2  Displacement  of  the  Free  End  of  the  Rod 

In  this  final  section,  we  find  expressions  for  the  six  scalar  displacements  (three 
translations  and  three  rotations)  of  the  &ee  end  of  the  rod  using  Castigliano’s  Second 
Theorem  which  may  be  stated  briefly  as:  If  a  linearly  elastic  structure  is  subjected  to  a  set 
of  loads,  then  the  displacement  of  a  load  in  its  direction  is  equal  to  the  derivative  of  the 
complementary  strain  energy  of  deformation  of  the  structure  with  respect  to  that  locui.  We 
may  write  this  mathematically  as: 

(3.10) 

Differentiating  the  complementary  strain  energy  expression  of  Appendix  A  with 
respect  to  Fx,  Fy,  F^,  Mx^  re^P^tively,  we  obtain  the  three  end  translations 

A  X,  A  K,  A  Z,  and  the  three  end  rotations  ^9x ,  Atfy ,  A^^  as: 
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—  hi^x^ fu^x'^ fn^Y  +  fit^z 

—  /21  ■'■  /22  -^y  +  /js  +  /24  ^x  +  /25  ^Y  +  fit  ^z 

=  fz\Fx^  hi^Y^  hz^z"^  h<  ^x  +  /s5  ^y  +  fzt  ^z 

^^x  =  /41  Fx^  UiFy^  /az^z-^  /aa  ^x  +  As  ^y  +  Ae  ^z 
=  hi  Fx^  hiFy-^  hz  Fz  +  fzA  ^x  '^hz^Y'^  ht  ^z 

—  hi  Fx'^  ftiFy'^  ftzFz"^  ftA  ^x  ■*■  hz  ^y  +  ht  ^z  (S-H) 

where  /^  represent  the  elements  of  the  6^6  symmetric  flexibility  matrix.  The 
complete  expressions  for  these  elements  are  listed  below.  For  the  sake  of  conciseness,  the 
following  notation  is  used:  c  =  cos<^,  3  =  sin  C  =  a/L,  and  S  =  b/L. 


{  6C*[Zac- ia^2(t>3^  +  <i>]  +  S*[Z<f>  +  2(f>^-Zsc] 

12  El  \  0 

-  6C*sM43-  33c-<^i  +  S*[-  3<^  +  2(^’  +  33c]  }  ^ 

r  3  /  ^2  c2  \ 

-  — —  i-15<i>-124>a^-2<l>^^A8a-33ac] 

UGJp\4>  J 

-  ——  (  9c\\  +  — ^  r  —  [(^  +  3C I  -r  2  5*  V3.I2) 

2EA  \<t>  ‘  /  AGA  \  0  ^  )  ’ 


—  [  {  2C*[Za'^-20ac  +  2c-2]^S^[0^- a^\ 

A E 1  \  0^ 

+  2C*S*[-  Za^  +  20a- 2c  +  2j +  S*[s*-  <^*]  }  ^ 
“  [8“  8c- ll3*  +  44>8c  +  4«^5- ^*]  ) 

2EA  \  0  ^  V  4Gi4  V  ^  / 


(  C*5* , 

8G7-  V  0^  ‘ 


(3.13) 
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/l3 


fu 


ho 


he 


hi. 


8GJ^ 


(3.14) 


(  ^  {  CS[4>a-4>^c]  +  2C^S[-2  +  2c-  2s^-^i4>3] 

\<i> 

-p  CS^  [4~  4c~  ^s~  tl>^c]  }  ^ 

(  [<i>^c  +  7</>s- 4a*  +  8c-  8]  ^ 

(  — 2  {  5[<6~3Cj  +  2C^5[^-)-5c~25)  +  5^[ac~<^!  } 

A  E I  \  (p  '  J 

■  8^7,  (^'“■'*’"“^0  (“> 

'•’  (M6) 


8GJ„ 


—  c  0»j-r  C5*[c- 1]  }  j 


12£?/ 


(  — 5  {  6C?^  [<!>  3ac  +  24>c^  ]  +  6C*5*  \  4>  +  Zac~  4^cl 
\  <t> 

■t"  5  ^  [  ~  3^  +  2^  *  +  Zac  ]  +  5*[3^  +  2^®~  Sw  ]  }  ^ 
3  /  o3  \ 

—  (-9<^-2^’-2Ve  +  l2^«*  +  33w]  j 


24  GJ. 


19 


fis  - 


4£/ 


+ 


(  {  C5:^c-  $- <f>^s]  +  2C^S[s  +  29c-  Z4>c\ 

V  4> 

-r  CS^[-  ZS- <f>C- <i>^9+i4>]  }  ^ 

EA\<t>  )  2GA.\  <f>  j 


(3.19) 


/24  - 

4E1 

{  S[-  a 

t-2c-  1 

] 

- 

53[,2. 

-*’i)) 

- 

8GJp 

/  C'^S 
V  <(>^ 

[  -  4-^4ci-a*'*-^*j 

J 

) 

(3.20) 

/j5  = 

L* 

4E1 

{  5[ac- 

'  <p]  2C^S  ~  scji- 

■5*I<^- 

sc]  } 

) 

- 

8GJp 

/  G’5 
V 

[^-  af]  ^ 

(3.21) 

he  - 

I*  1 
El  ' 

C^\<t>c 

-aj  +  C5’ja-.^;  }  ^ 

1 

- 

I* 

2GJp 

f  CS^ 

\ 

[-  ^c-  <^  +  2a]  ^ 

(3.22) 

fiz  = 

2E1 

{  C^[4> 

—  ac  ]  +  G  *  ( 3^  +  St 

•~  4#]  } 

) 

-  f 

4GJp  V 

c* 

1 

1 

«• 

o 

IS’l  + 

L 

2GA 

[C*K3-23) 
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/«  =  ci*.i+csMJ«-J+^»l}) 

/ss  ~  YWl  C7[5-^c!  +  C5*(a-^c]}  ^ 


fae  - 

f  C 
El\<p 

*5.  I*  / 

i«-^i  ) 

(3.26) 

/44  = 

i  (1 

2EI  \ 

^  ( 
iGJ  V 

(3.27) 

/4S  = 

{Ar 

(3.28) 

fit  = 

(-- 
V  El 

(3.29) 

/ss  = 

^  (' 
2EI  \  ' 

^  1 
AGJ  ' 

(3.30) 

ht  = 

(  " 
\2GJ 

(3.31) 

/ee  = 

-t  IC» 
El  ‘ 

1+  ^  [5»1 

'  2GJ‘  ' 

(3.32) 

Equations  (3.11)  represent  the  scalar  displacements  and  rotations  of  the  end  of  the 
helical  rod,  in  terms  of  the  XYZ  components  of  the  applied  force  and  moment,  and  the 
geometrical  parameters  defining  the  helix,  namely,  a  and  By  making  suitable 
substitutions,  it  is  possible  to  reduce  them  to  those  of  a  cantilevered  straight  beam.  This  is 
demonstrated  in  the  following  subsection. 
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3.2.1  Reduction  of  the  General  Equations  to  the  Straight  Beam  Case 


The  equations  for  the  helical  arc  may  be  reduced  to  those  for  the  straight  beam  in 
two  equivalent  ways;  by  setting  a  =  0,  which  forces  ^  to  be  zero,  and  alternatively  by 
letting  <t>  approach  zero.  This  is  illustrated  in  Figure  3.1.  It  is  important  to  realize  that  we 
cannot  simply  substitute  ^  =  0  into  the  general  equations,  but  must  express  sin0  and  cos^ 
as  Taylor  series  about  ^  =  0  and  then  neglect  the  higher  order  terms.  Both  the  methods 
yield  identical  results,  which  are  exactly  the  same  as  the  strength-of-materials  deflections 
for  the  straight  cantilevered  beam.  We  present  below  the  reduced  equations  that  result 
from  applying  the  first  method: 


Air  = 

FxL^  FxL 

ZEI  2GA 

(3.33) 

c> 

II 

FyL^  FyL 

ZEI  2GA 

(3.34) 

A^  = 

FzL 

EA 

(3.35) 

= 

FyL^ 

2EI 

(3.36) 

A0y  = 

FxL^ 

2EI 

(3.37) 

— 


0 


(3.38) 
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S.2.2  Reduction  to  a  Circular  Arc  of  Small  Curvature  in  the  XY  Plane 

The  reduction  of  the  general  equations  to  the  case  of  a  circular  arc  in  the  XY 
plane  is  accomplished  by  setting  6  =  0  or  a  =  £,  and  expressing  sin^  and  cos^  as  Taylor 
series  about  ^  =  0  2md  retting  terms  upto  the  order  of  <f>^.  We  have,  then,  the  following 
results: 


Note  that  if  we  now  set-to  zero  in  the  above  equations,  we  obtain  the 
displacements  of  a  straight  cantilever  beam  in  the  XY  plane,  analogous  to  those  presented 
for  the  case  in  which  a  =  0,  presented  in  Section  3.2.1.  In  the  next  chapter,  we  impose  two 
different  conditions  of  compatibility  of  displacements  of  the  muzzle  ends  of  the  gun  tube 
and  the  shroud.  We  thus  obtain  the  force  exerted  on  the  gun  tube  by  the  thermal  shroud 
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due  to  its  expeuision,  and  thereby,  the  displacement  of  the  end  of  the  gun  tube. 


CHAPTER  4 


IMPOSITION  OF  COMPATIBILITY  CONSTRAINTS 

4.1  A  Convenient  Coordinate  System 

In  chapters  2  and  3,  we  developed  an  analytical  solution  for  the  displacement  at 
the  end  of  a  helical  rod  of  constant  circular  cross-section  in  terms  of  an  applied  force  and 
moment.  However,  in  the  actual  problem,  we  do  not  know,  a  priori,  what  the  value  of  the 
force  exerted  on  the  gun  tube  by  the  shroud  is.  The  problem  is  statically  indeterminate. 
Hence  we  must  rely  on  the  imposition  of  compatibility  of  the  displacements  of  the  ends  of 
the  gun  tube  and  the  thermal  shroud  to  resolve  this  indeterminacy. 


Z 


Figure  4-1:  The  Chord  System 
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Note  that  by  varying  the  parameters  a  and  (f>,  we  can  vary  the  curvature  and 
torsion  of  the  helical  rod.  In  order  to  perform  a  parametric  study  of  the  displacement,  we 
must  vary  the  curvature  k  and  torsion  A  of  the  helical  rod,  and  solve  the  compatibility 
conditions  that  we  impose  for  each  such  case  to  obtain  the  force  acting  on  the  rod,  and  then 
compute  the  six  displacements  of  the  muzzle  end.  For  this  purpose,  we  make  the 
reasonable  assumption  that  the  centroidal  axis  of  the  thermal  shroud  passes  through  the 
two  end  points  A  and  B  of  the  helical  rod,  which  define  a  chord  between  them.  This  is 
shown  in  Figure  4.1.  We  therefore  define  a  coordinate  system  whose  Y  axis  is  in  the 
direction  of  the  vector  AB,  and  whose  X  axis  lies  in  a  plane  parallel  to  the  XY  plane,  so 
that  if  we  set  6  =  0  and  let  o  approach  zero,  the  coordinate  system  becomes  identical  to  the 
global  coordinate  system.  This  coordinate  system  is  called  the  chord  system. 

The  chord  system  is  the  coordinate  system  that  we  shall  use  in  Chapter  5  for  the 
parametric  studies  that  compare  displacements  for  different  initial  curvatures  and  torsions 
of  the  gun  tube,  and  different  axial  and  bending  sti&esses  of  the  shroud. 

4.2  Unit  Vectors  in  the  Chord  System 

Let  X^Y^Z^  denote  the  chord  system  shown  in  Figure  4.1.  The  position  vectors  of 
points  A  and  B  are: 

^{A)  =  i  i 

<P 

p{B)  =  -  COS0  i  -  sin0  j  +  6k  (4.1) 

4> 

Therefore,  the  vector  AB  is  given  by: 

AB  =  -  (cos^ -  l)l  +  -mn^j  +  6k 
<j>  (j) 

From  Equation  (4.2),  we  may  write  the  unit  vector  along  the  Y^  axis  as: 


(4.2) 
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I  _  AB  _  (cos<t>  -  l)i  +  sinc^j  +  (06/a)k 
‘  '■  AB  V  2(l-cos^)  -r  {<pb/a)^ 


(4.3) 


Having  obtained  j  in  terms  of  i,  j  cmd  k,  the  unit  vectors  in  the  global  coordinate 

C 

system,  we  now  proceed  to  obtain  and  k^.  As  mentioned  in  Section  4.1,  the  axis  is 
parallel  to  the  glob£d  XY  plane.  Thus  may  be  written  as: 


i  =  pi  -  ^  Ok 


(4.4) 


Since  and  are  mutually  perpendicular,  their  scalar  product  must  be  zero,  i.e.. 


i  .j  =  0 

c  "c 


(4.5) 


That  is. 


p  ( cos  0  -  1 )  —  5  sin  <?  =  0.  (4.6) 

Since  i  is  a  unit  vector,  we  have: 
c 

=  l  (4.7) 

Solving  equations  (4.6)  and  (4.7)  for  p  and  we  get 

i  =  i  -^  ( 1  ~  cos<p  )  j 
V  2  ( 1  -  cos  0  ) 

Having  obtained  1^  and  j^,  we  may  simply  write  k^  as; 

lj  =  i  X  j  =  (</>6/a)  ( 1-  co8<^)  i  -  (i^6/a)  sin(^j  +  2(1-  co8<^)k 
'  '  ‘  \/4(l-  cos^)*  -t  2(^6/o)*(l  -  cos^) 

Equations  (4.3),  (4.8)  and  (4.9)  give  us  the  desired  unit  vectors  in  the  chord 
system  in  terms  of  the  unit  vectors  in  the  global  XYZ  system.  From  these  three  equations, 
we  now  derive  the  chord-to-global  transformation  matrix,  CG,  that  relates  the  components 
of  displacement  and  force  in  the  chord  system  to  their  equivalents  in  the  global  system. 
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4.8  CG:  The  Chord-to-Global  Transformation  Matrix 


Consider  a  generic  vector,  u^,  that  is  to  be  transformed  from  the  chord  coordinate 

system  into  the  global  coordinate  system.  Let  s  be  the  matrix  whose  rows  consist  of  the 

coefficients  of  i,  j  and  k  unit  vectors  in  equations  (4.3),  (4.8)  and  (4.9)  respectively.  If  is 

written  as  u  =  u.,  i  *  u.,,  j  ^  u_  k  ,  then  substituting  for  i  ,  J  and  k  from  the  above- 
c  Xc  c  Yc  c  Zc  c  °  c  c  c 

mentioned  equations,  we  get 


s 

i 

3  _  j  3..  k 

Xc 

11 

12*'  13 

ti,. 

,5 

i  — 

3.  j  3  .  k 

Yc 

21 

22*'  23 

i  — 

3.  j  -r  3„  k 

Zc  • 

31 

32*'  33 

Rearranging  terms,  we  get  the  following: 


(4.10) 


u 

c 


^12  “Xe  ~  *22  *32  ’‘Zc  '  ^ 


*13  ^Xc  ~  *23  “Yc  *33  “zc  '  ^ 


(4.11) 


Clearly,  from  Equation  (4.11),  we  see  that  CG,  the  chord-to-global  transformation  matrix 
is  the  tr2uispose  of  s: 


CG  =  8^ 


(4.12) 


f  T 

Thus,  to  convert  the  components  of  a  vector  u  =  ;  ]  from  the  chord 

system  to  the  global  system,  we  pre-multiply  u  with  CG.  The  global-to-chord 

transformation  is  achieved  by  inverting  CG.  But  since  the  transformation  is  orthogonal, 
“1  7* 

GC  =  CG'  =  CG  =  8.  In  Section  2.1,  we  presented  the  expressions  for  the  tangent, 
binormal  and  principal  normal  to  a  helix,  at  a  generic  point  along  its  length,  in  terms  of  the 
global  unit  vectors  i,  j  and  k.  The  matrix  formed  by  the  coefBcients  of  these  tmit  vectors  in 
equations  (2.7),  (2.8)  and  (2.9)  similarly  represent  the  elements  of  the  global-to-tangent 
transformation  matrix  at  any  pioint  along  the  helix. 
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4.4  Flexibility  Matrix  in  the  Chord  System 

Since  the  position  of  the  thermal  shroud  in  the  chord  system  is  always  fixed  while 
it  varies  in  the  global  coordinate  system  with  changing  curvature  and  torsion,  it  is  clear 
that  compatibility  of  the  displacements  of  the  ends  of  the  gun  tube  and  the  shroud  must  be 
imposed  in  the  chord  system.  This  means  that  we  must  transform  the  flexibility  matrix  of 
Equation  (3.11),  which  yields  displacements  in  the  global  system,  to  one  that  gives 
displacements  directly  in  the  chord  system. 

An  important  point  to  note  here  is  the  assumption  that  the  £Lngular  displacements 
and  may  be  decomposed  into  components  along  the  coordinate  axes  emd 
transformed  between  different  coordinate  systems  in  exactly  the  manner  the  translational 
displacements  are.  This  is  acceptable  for  small  angular  displacements  as  in  our  case. 


Consider  a  6  x  6  chord-to-global  tremsformation  matrix  T,  which  is  block 

diagonal,  both  the  3  x  3  diagonal  blocks  being  equal  to  CG.  Let  f  and  f  denote  the  6x6 

flexibility  matrices  in  the  global  and  chord  coordinate  systems.  Further,  let  A,  F  and  A  , 

c 

F^  represent  the  6x1  displacement  and  force  vectors  in  the  global  and  chord  systems 
respectively.  Then  we  have  the  following: 


[> 

II 

(4.13) 

^  =  f  F 

c  c  c 

(4.14) 

A  =  TA 

c 

(4.15) 

f  =  Tf 

c 

(4.16) 

Substituting  equations  (4.15)  and  (4.16)  into  Elquation  (4.13),  and  pre-multiplying  both 
sides  by  T”^  =  T^,  we  get: 

A  =  'T^fTiF 

C  •  ‘  C 


(4.17) 
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Comparing  equations  (4.14)  and  (4.17),  we  get  f  as: 

f  =  T^fT  (4.18) 

4.5  Compatibility  1:  Shroud  with  Zero  Transverse  StifEaess 

We  assume  a  condition  of  compatibility  of  displacements  of  the  thermal  shroud 
and  gun  tube  ends  wherein  the  shroud  has  stifhiess  only  in  its  axial  direction.  That  is,  the 
shroud  has  no  transverse  or  rotationed  stiffiiess.  Thus  translational  movement  is  permitted 
to  the  gun  tube  end  freely  in  the  transverse  direction,  perpendicular  to  the  centroidal  axis 
of  the  shroud.  As  a  result,  there  is  no  opposing  force  exerted  by  the  shroud  in  the 
transverse  direction.  In  fact,  the  only  force  exerted  by  the  shroud  on  the  gun  tube  acts 
along  the  axis  of  the  chord  system.  Rotations  of  the  gun  tube  end  are  also  permitted 
about  all  three  axes  in  the  chord  system,  and  similarly,  there  are  no  moments  exerted  at  the 
gun  tube  end. 

If  is  the  force  exerted  on  the  end  of  the  gun  tube,  then  its  displacement  in  the 
direction  of  the  force  is  given  as 

Ar(^nft*6«)  =  (4-19) 

where  /^22  ^  element  belonging  to  the  2”*^  row  and  2**'^  colunon  of  the  chord  flexibility 
matrix.  The  gun  tube,  by  Newton’s  Third  Law,  must  in  its  turn  exert  an  equal  but 
opposite  force  on  the  end  of  the  shroud.  Let  a.  Eg,  A,  and  A  7  be  the  length,  coefficient 
of  thermal  expansion,  modulus  of  elasticity,  cross-sectional  area  and  rise  in  temperature  of 
the  shroud  respectively.  The  displacement  of  the  end  of  the  shroud  is  composed  of  two 
contributions,  namely,  the  thermal  expansion  and  the  contraction  due  to  the  opposing  axial 
force  exerted  by  the  gun  tube.  Thus  we  have 

Fy^L. 

A  Y  {shroud)  =  A  7  -  - 

*  EgAg 


(4.20) 
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Compatibility  asserts  that  the  displacements  in  equations  (4.19)  and  (4.20)  be  equal.  Thus 
we  find  the  force  acting  on  the  end  of  the  gun  tube  to  be 


F  = 

''' 

where  L,,  the  length  of  the  shroud,  may  be  easily  seen  from  Figure  2.1  to  be 


(4.21) 


L,  =  i  '  -  ( cos  <?>  -  1 )  ‘  *  -(-  !  -  sin  6  *  ^  1 .  (4.22) 

The  six  displacements  of  the  end  of  the  gun  tube  may.  now  be  found  by  substituting  f 
into  equations  (3.11)  rewritten  with  elements  of  the  flexibility  matrix  in  the  chord  system. 


4.6  Compatibility  2:  Shroud  with  Finite  Non-Zero  Transverse  StifBaess 


In  this  section,  we  consider  a  different  kind  of  compatibility  condition  in  which  we 
assume  the  shroud  to  have  a  finite  non-zero  transverse  sti&ess.  The  gun  tube  end  is 
permitted  both  transverse  displacements  as  well  as  rotations  relative  to  the  shroud. 
However,  since  the  shroud  now  has  sti&ess  in  the  transverse  direction,  it  will  exert  a 
resistance  to  the  transverse  displacement  of  the  gun  tube  end.  An  equivalent  way  of 
representing  such  a  compatibility  relationship  is  to  consider  the  gun  tube  end  as  being 
attached  to  the  shroud  end  through  a  ball  and  socket  joint. 

Since  the  shroud  now  offers  resistamce  to  the  transverse  motion  of  the  gun  tube 
end,  f  and  F are  non-zero.  The  compatibility  of  all  three  components  of  displacement 
of  the  ends  of  the  gun  tube  and  the  shroud  must  be  enforced  in  order  to  determine  the 
forces  on  the  gun  tube.  Thus,  the  imposition  of  compatibility  results  in  a  set  of  three 
simultaneous  equations  iaF,F  and  F„ ,  which  are  presented  below: 

Xc  ic  Zc 
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Fx^W 

~  TeJ,  "  ^Yc  ^ 

^Yc^> 

"« ■^# 


^.cV 


-  -^cSl  ^Jfc  ^  4»2^yc  ^  /c»3^Ze 


(4.23) 


Here,  /,  represents  the  moment  of  inertia  of  the  cross-section  of  the  shroud.  It  is  worth 
noting  that  if  the  inertia.  of  the  shroud  were  to  approach  zero  (which  would,  in  physical 
terms,  amount  to  the  shroud  becoming  infinitely  flexible),  F  and  F  ,  which  represent  the 
opposition  of  the  shroud  to  bending,  would  also  tend  to  zero,  and  the  above  compatibility 
conditions  would  reduce  to  those  of  the  first  case.  The  next  chapter  contains  peu^ametric 
plots  showing  the  variation  of  the  forces  and  displacements  experienced  by  the  end  of  the 
gun  tube  for  different  initieil  curvatures  and  torsions,  for  both  of  the  above  compatibility 
conditions. 


CHAPTER  6 


A  PARAMETRIC  STUDY  OF  DISPLACEMENT 

6.1  Introduction 

Having  obtained  a  closed-form  solution  for  the  end  displacement  of  the  gun  tube 
in  Chapter  3,  <ind  developed  equations  of  compatibility  of  gun  tube  and  shroud 
displacements  in  Chapter  4,  we  proceed  to  perform  a  par£tmetric  study  of  the  force  and 
displacement  experienced  by  the  end  of  gun  tube.  The  parameters  in  this  study  ue  the 
initial  curvature  and  torsion  of  the  gun  tube,  the  kind  of  compatibility  condition  imposed, 
namely,  one  of  the  two  kinds  discussed  in  sections  4.5  and  4.6,  and  the  axial  and  flexural 
rigidity  of  the  thermal  shroud. 

Several  plots  of  force  and  displacement  components  m  the  chord  system  are 
included  in  this  chapter,  usually  accompanied  by  discussions  regardmg  their  physical 
interpretation.  These  are  ctdculated  for  a  fixed  rise  in  temperature  of  100  ‘  C  of  the  thermal 
shroud.  The  variation  of  the  components  of  force  and  displacement  with  temperature 
would,  of  course,  be  linear.  The  gun  tube  and  shroud  data  assumed  in  the  computations 
follow. 

Gun  Tube: 

Length  =  5.25  m 

Inner  radius  =  0.06  m 

Outer  radius  =  0.10  m 

Elasticity  modulus  =  4.35097  x  10*  N/m* 

Shear  modulus  =  1.67345  x  10*  N/m* 
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Thermal  Shroud: 

Inner  radius  =  0.11  m 

Outer  radius  =  0.145  m 

Elasticity  modulus  =  4.35097  x  10®  N/m* 

Coefficient  of  thermal  expansion  =  1.206  x  10~®  /  *  C 

Each  parametric  plot  is  a  surface  depicting  the  variation  of  the  particular  force  or 
displacement  component,  taken  in  the  chord  system,  for  different  values  of  initial  curvature 
and  torsion  of  the  gun  tube.  It  is  assumed  that  the  length  of  the  gun  tube  is  a  constant. 
The  length  of  the  thermal  shroud  is  the  distance  between  the  two  ends  of  the  gun  tube,  and 
therefore  varies  with  the  gun  tube  geometry  as  given  by  Elquation  (4.22).  The  force  and 
displacement  components  are  plotted  for  the  normal  values  of  axial  rigidity,  EA,  and 
flexureil  rigidity,  El,  of  the  thermal  shroud,  and  also  for  high  and  low  EA  and  El,  the  values 
being  increased  and  decreased  respectively  by  a  fturtor  of  a  thousand.  In  the  following 
section,  we  define  the  geometrical  parameters  employed  in  the  study,  namely,  the  initial 
curvature  and  torsion  of  the  gun  tube,  in  terms  of  the  corresponding  geometrical 
parameters  that  define  the  shape  of  the  helix  of  Figure  2.1,  namely,  a  and  4>. 

6.2  The  Definition  of  Ctirvature  and  Torsion 

Every  space  curve  of  constant  curvature  and  torsion  is  uniquely  and  completely 
specified  by  these  two  quantities,  apart  from  its  position  in  space.  While  it  would  be 
possible  to  consider  a  and  <(>,  which  define  the  shape  of  the  helix  of  Figure  2.1,  as  a 
legitimate  choice  of  geometrical  parameters  for  this  study,  we  have  chosen  to  employ  the 
initial  curvature,  k,  and  torsion,  A,  of  the  centroidal  axis  of  the  gun  tube  as  parameters  to 
represent  its  initial  shape  on  account  of  their  generality.  We  now  derive  relationships 
between  these  two  sets  of  geometrical  parameters. 
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The  curvature  2ind  torsion  of  a  space  curve  are  defined  by  the  following  two 
Frenet-Serret  formulae, 


—  =  Kl/,  K  >  0  (5.1) 

ds 


da 


(5.2) 


where  r,  0  and  v  are  the  vector  tangent,  binortnal  and  principal  normal  to  the  curve,  and  a 
is  the  eirc  length.  By  definition,  k  and  A  are  zero  for  a  straight  line.  The  cross  product 
r  X  V  defines  the  binormal  vector,  0,  which  is  the  outward  normal  to  the  osculating  pltme. 
The  torsion.  A,  measures  the  rate  at  which  the  direction  of  0  changes  along  the  curve. 
From  equations  (2.7),  (2.8),  (2.9),  (5.1)  and  (5.2),  the  curvature  and  torsion  of  the  helix  are 
found  to  be: 


K 


(5.3) 


A 


6b  ab 


(5.4) 


The  parameters  a,  6,  R  and  L  in  equations  (5.3)  and  (5.4)  are  exactly  as  shown  in  Figure 
2.1.  Note  that  R  represents  the  radius  of  the  cylindrical  surface  that  contains  the  helix. 

Note  also  that  the  length  of  the  gun  tube,  L,  is  fixed,  and  6  is  given  as  V  -  o* . 


Since  we  are  only  concerned  with  the  magnitude  of  torsion,  for  the  purposes  of  this 
analysis,  we  rewrite  Equation  (5.4)  as 


A 


From  equations  (5.3)  and  (5.5),  we  get 


(5.5) 


K  _  a 

X  ~  ~b’ 


(5.6) 


and  thus,  we  have  the  following  two  expressions  for  a  and  b  in  terms  of  k  and  A: 
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a 


kL 


V  K  *  -  A  * 


(5.7) 


h  = 


(5.8) 


From  equations  (5.3)  and  (5.7),  we  also  obtain  the  following  expressions  for  ^  and  R  in 
terms  of  k  and  A: 


<t>  =  L  V  K  *  —  A  * 


(5.9) 


R  = 


(5.10) 


Having  obtained  the  curvature,  k,  and  torsion,  A,  in  terms  of  a,  tp  tmd  constant  L, 
and  also  developed  reciprocal  relationships  for  a,  i,  <P  and  R  in  terms  of  k,  A  and  L,  we  are 
equipped  to  compute  the  force  and  displacement  components  that  we  seek  to  plot,  given 
the  values  of  initial  curvature  and  torsion  of  the  gun  tube.  This  may  be  done  by  first 
computing  a,  b,  tp  and  R  from  equations  (5.7),  (5.8),  (5.9)  and  (5.10)  resp>ectively,  and 
from  Ek^uation  (4.22),  then  solving  the  particular  compatibility  equations,  equations  (4.21) 
or  (4.23)  as  the  case  may  be,  to  obtain  the  applied  force  components,  and  lastly, 
substituting  these  scalar  forces  into  equations  (3.11)  to  get  the  required  displacements  and 
rotations  of  the  gun  tube  end. 


5.8  Compatibility  1:  Physical  Interpretation 

In  Compatibility  1,  we  assume  that  the  shroud  has  sti&ess  only  in  its  axial 
direction.  The  shroud  is  assumed  to  have  zero  transverse  or  rotational  sti&ess.  Since  the 
gun  tube  end  is  free  to  move  in  the  transverse  direction,  the  only  force  that  acts  on  it  is 
directed  along  the  centroidal  axis  of  the  shroud. 
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6.3.1  ’’Iii'Plane’’  Force  and  Displacements:  An  Anomaly  Explained 

Each  point  on  the  curve  has  an  associated  osculating  plane  which  is  normal  to  the 
binormal  vector  at  that  point.  The  variation  of  the  osculating  plane  along  the  curve  for  the 
values  of  curvature  and  torsion  considered  is  very  small.  In  other  words,  the  helical  arc 
described  by  the  centroidal  axis  of  the  gun  tube,  and  the  chord  along  which  the  shroud  is 
2digned.  form  a  surface  that  is  almost  planar.  The  non-planarity  of  this  surface  b  measured 
in  terms  of  the  rate  of  change  of  the  binormal  vector  along  the  arc,  A  =  |  djd/ds  | ,  and  is  at  a 
maximum  when  the  curvature  and  torsion  are  equal.  This  is  illustrated  in  Figure  5.1,  for  a 
fixed  radius.  R.  It  can  be  proven  mathematically  that  this  characteristic  is  unchanged  even 
if  the  radius  is  allowed  to  vary  along  a  constant  curvature  contour. 

This  surface  is  almost  identical  to  the  plane  in  the  chord  system  of  Figure 
4.1.  The  forces  and  displacements  in  the  X^Y^  plane,  namely,  Fj^g,  Fy^,  and  V^, 
and  62c:  may  be  regarded  as  being  in  the  plane  of  the  curve  and  will  be  referred  to  in  the 
following  discussion  as  the  ’’in-plane”  forces  and  displacements.  Since  all  the  forces  and 
displacements  we  refer  to  are  in  the  chord  system,  we  will  drop  the  "chord”  subscript  for 
convenience. 

Fx  and  F2  are  identically  zero  due  to  the  fact  that  the  shroud  offers  no  resistance 
to  the  motion  of  the  gun  tube  end  in  the  transverse  direction.  Figure  5.2  shows  a  plot  of 
Fy.  Note  that  Fy  decreases  as  the  curvature  increases,  whereas  it  is  virtually  unaffected  by 
torsion.  The  decrease  of  Fy  with  increase  in  curvature  is  consistent  with  experience  with 
rods  with  plane  curvature.  The  apparent  independence  of  Fy  on  torsion  is  somewhat 
surprising  at  first  sight,  since  both  curvature  and  torsion  may  be  thought  of  as  representing 
deviations  from  straightness.  At  a  physical  level,  this  may  be  rationalized  by  saying  that 
the  helical  arc  representing  the  centroidal  curve  of  the  gun  tube  is  an  almost  circular  arc  in 
the  X^Y^  plane  whose  curvature  b  held  fixed  (along  a  constant  curvature  contour).  Thus 
increasing  torsion  only  mcreases  the  deviation  of  the  curve  from  the  X^Y^  plane  upto  a 
maximum  at  k  =  A.  As  A  grows  greater  than  k,  the  curve  approaches  a  straight  line,  since 
h  becomes  much  greater  than  a.  Thb  leads  us  to  expect  an  increase  in  Fy  with  increasing 
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torsion.  This  is  borne  out  by  the  numbers  generated  for  the  plot  of  Figure  5.2.  The 
increase  in  Fy  is  more  apparent  for  low  values  of  curvature  than  for  high  values,  because 
when  the  curvature  is  small,  Fy  is  very  close  to  the  limiting  value  for  a  straight  rod. 

The  behavior  observed  for  Fy  also  occurs  in  the  plots  of  U,  V  and  0^  (figures  5.3  • 
5.5),  namely,  they  do  not  show  an  appreciable  variation  with  torsion  along  lines  of  constant 
curvature.  The  same  argument  holds,  namely,  the  torsion  changes  the  deviation  of  the 
helical  arc  from  the  plane  but  the  curve  is  essentially  circular,  has  the  same  curvature 
and  is  almost  contained  in  the  X  Y  plane,  so  that  the  displacements  are  virtually 

C  C 

unaffected  by  a  variation  of  torsion.  Note  that  is  always  negative;  this  demonstrates  that 
the  helical  rod  “unwinds”  as  it  is  stretched. 


The  U  and  0^  displacements  show  a  maximum  along  lines  of  constant  torsion. 
This  is  a  result  of  the  imposition  of  compatibility  and  may  be  explained  as  follows:  Initifdly, 
the  curvature  lends  the  rod  an  added  flexibility  as  a  result  of  which  the  displacements 
increase.  However,  as  curvature  increases,  the  force  developed  due  to  the  interaction 
between  the  shroud  and  the  gun  tube,  Fy,  decreases  monotonically.  A  point  is  reached 
when  the  decrease  in  the  force  overtakes  the  increase  in  flexibility. 


5.3.2  Out-of-Plane  Displacements 


The  W,  Ox  ^d  Oy  displacements  of  figures  5.6  -  5.7  {Ox  is  negligibly  small  and  has 

not  been  shown),  on  the  other  hand,  represent  motions  of  the  end  of  the  gun  tube  in  a 

direction  perpendicular  to  the  X  Y  plane.  These  are  therefore  termed  as  the  '’out-of- 

c  c 

plane”  displacements.  Ox  remuns  essentially  zero  for  all  values  of  the  applied  force.  Since 

the  applied  force  is  in  the  X  Y  plane,  and  the  helical  arc  is  also  essentially  confined  to  this 

c  c 

plane,  we  naturally  expect  that  the  out-of-plane  displacements  will  be  smaller  in 
magnitude.  This  is  seen  to  be  the  case.  These  displacements  show  peaks  in  the  region 
defined  by  k  =  A.  This  is  because  when  x  =  A,  the  rod  departs  from  the  plane  the  most.  In 
the  plot  of  W  however,  the  peak  is  off  the  k  =  A  line  of  symmetry.  This  may  be  attributed 
to  the  imposition  of  compatibility  as  explained  in  Section  5.2.1. 
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AXIS 

X  (Curvature) 

Y  (Torsion) 

Z  (Dependent  Venable 
plotted  for  a  temp, 
rise  of  100  deg.  C) 


SCALE 

1  unit  -  Cur  X  10.0  (Log) 
1  unit  -  Tor  X  10.0  (Log) 
1  unit  -  13440.0000  N 


RANGE 

l.OE-4  /■  to  1,0E-1  /■ 

0.  l,0E-4  /n  to  l.OE-1  /m 
5190.7310  N  to  61442.1406  N 


FOflCEY 


TOR 


CUR 

Compatibility  1:  Fy  (Chord)  vs  Cur  vs  Tor 


Figure  6-2:  Compatibility  1:  Fy  versus  Curvature  versus  Torsion 
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AXIS 

X  (Curvature] 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  C) 


SCALE 

1  unit  >  Cur  X  10.0  (Log) 
1  unit  •  Tor  x  10.0  (Log) 
1  unit  •  0.0420  ffl 


RANGE 

l.OE-4  /«  to  l.OE-1  /» 

0,  l.OE-4  /«  to  l.OE-1  /» 
0.0006  m  to  0.0990  a 


U 


Compatibility  1:  U  vs  Cur  vs  Tor 


Figure  6-S:  Compatibility  1:  U  versus  Curvature  versus  Torsion 
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AXIS 

X  (Curvature] 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  C) 


SCALE 

1  unit  ■  Cur  X  10.0  (Log) 
1  unit  «  Tor  x  10. 0  (Log) 
1  unit  -  0.0024  IB 


RANGE 

l.OE-4  /IB  to  l.OE-1  /« 

0.  l,0E-4  /fli  to  l.OE-i  /m 
0.0037  m  to  0.0060  m 


Compatibility  1:  V  vs  Cur  vs  Tor 


Figure  6-4:  Compatibility  1:  V  versus  Curvature  versus  Torsion 
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AXIS  SCALE 

X  (Curvature)  1  unit  •  Cur  x  10.0  (Log) 

Y  (Torsion)  1  unit  -  Tor  x  10.0  (Log) 

Z  (Dependent  Variable  l  unit  0.4800  deg 

plotted  for  a  teop. 
rise  of  100  deg.  C) 


RANGE 

l.OE-4  /■  to  l.OE-1  /■ 

0,  l.OE-4  /•  to  l.OE-1  /■ 
-2. 1629  deg  to  -0.0142  deg 


Compatibility  1:  Thetaz  vs  Cur  vs  Tor 


Figure  B-6:  Compatibility  1;  versus  Curvature  versus  Torsion 
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Note  that  in  Compatibility  1,  the  displacement  is  dominated  by  U,  that  is,  the  U 
component  of  displacement  is  of  the  largest  magnitude.  This  is  because  V  is  constrained 
due  to  the  sudal  rigidity  of  the  shroud,  and  W  is  small  because  it  is  an  out-of-plane 
displacement. 

Note  also  that  the  anomalous  behavior  with  respect  to  torsion  is  largely  due  to  the 
fact  that  the  displacement  components  plotted  are  in  the  chord  system.  In  other  words,  if 
the  XY  plane  of  the  system,  in  which  the  displacements  are  described,  were  not  to  almost 
contain  the  helical  curve,  there  would  no  longer  be  any  sense  in  which  a  displacement  could 
be  considered  to  be  in-plane  or  out-of-plane,  and  the  effects  of  these  two  kinds  of  motions 
would  be  fused. 

S.4  Compatibility  2  and  its  Variations 

In  Compatibility  1,  we  assumed  that  the  shroud  is  infinitely  flexible  in  the 
transverse  direction.  Compatibility  2,  on  the  other  htmd,  is  a  condition  wherein  the  shroud 
is  assumed  to  have  a  finite  non-zero  transverse  stii&ess.  This  gives  rise  to  two  additional 
components  of  force  exerted  on  the  gun  tube  at  its  end,  namely,  Fx  and  F^- 

6.4.1  Compatibility  2 

The  shape  of  the  plots  of  the  different  force  and  displacement  quantities  against 
curvature  and  torsion  is  unchanged  by  the  change  of  the  compatibility  condition.  In  other 
words,  the  in-plane  forces  and  displacements,  and  U,  V,  figures  5.8,  5.9  and 

5.11,  5.12,  5.16  respectively  show  the  same  lack  of  variation  with  torsion  along  lines  of 
constant  curvature.  In  addition,  U  (and  the  corresponding  force  of  reustance  Fx)  and  $2 
show  the  same  maxima  along  constant  torsion  contours.  As  explained  earlier,  these  are  due 
to  the  imposition  of  compatibility.  The  out-of-plane  displacements  W,  €xi  und  By  also 
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AXIS  SCALE 

X  (Curvature)  1  unit 

Y  (Torsion)  1  unit 

Z  (Dependent  Variable  1  unit 

plotted  for  a  teap. 
rise  of  lOO  deg.  C) 


-  Cur  X  10.0  (Log) 

■  Tor  X  10.0  (Log) 

■  0.0020  ■ 


RANGE 

S.OE-4  /m  to  l.OE-1  /• 

0.  l.OE-4  /■  to  l.OE-1  /■ 


0.0000  ■  to 


0.0017 


TOR 


\ 


\ 


\ 


CUR 


Compatibility  1:  W  vs  Cur  vs  Tor 


Figure  6-6;  Compatibility  1:  IV  versus  Curvature  versus  Torsion 


AXIS 

X  (Curvature) 

Y  (Torsion) 

Z  (Dependent  VariaDle 
plotted  for  a  tenp. 


SCALE  HANBE 

1  unit  -  Cur  X  10.0  iLog)  l.OE-4  /■  to  l.OE-1  /• 

1  unit  -  Tor  x  10. 0  (Log)  0.  l.OE-4  /■  to  l.OE-1  /■ 

1  unit  >  0.0030  deg  -0.0039  deg  to  0.0000  deg 

thetay 


Figure  B-T:  Compatibility  1: versus  Curvature  versus  Torsion 
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AXIS  SCALE 

X  (Curvature)  1  unit 

Y  (Torsion)  1  unit 

Z  (Dependent  Variable  i  unit 

plotted  for  a  temp, 
rise  of  100  deg.  C] 


RANGE 

Cur  X  10.0  (Log)  l.OE-4  /■  to  l.OE-1  /■ 

Tor  X  10.0  (Log)  0.  l.OE-4  /■  to  l.OE-1  /■ 

400.0000  N  -657.4645  N  to  -3.0925  N 

FORCEX 


Compatibility  2: 


Fx  vs  Cur  vs  Tor 


Figure  6-8:  Compatibility  2:  Fx  versus  Curvature  versus  Torsion 
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AXIS  SCALE  RANGE 

X  (Curvature)  1  unit  «•  Cur  x  10. 0  (Log)  i.0E**4  /■  to  i.OE-1  /■ 

Y  (Torsion)  1  unit  -  Tor  x  lO.o  (Log)  0.  l.OE-4  /■  to  l.OE-1  /■ 

Z  (Dependent  Variable  1  unit  -  13440.0000  N  9294. 3164  N  to  61442.4570  N 

plotted  for  a  temp, 
rise  of  100  deg.  C) 


FOflCEY 


CUB 

Compatibility  2:  Fy  vs  Cur  vs  Tor 


Figure  6-9:  Compatibility  2;  Fy  versus  Curvature  versus  Torsion 
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AXIS 

X  (Curvature] 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  tenp. 
rise  of  100  deg.  C) 


SCALE 

1  unit  •  Cur  X  10.0  (Log) 
1  unit  Tor  X  10.0  (Log) 
1  unit  «  11.0000  N 

F0RCE2 


RANGE 

l.OE-4  /m  to  l.OE-l  /II 
0.  l.Oe-4  /«  to  l.OE-l  /■ 
-11.4691  N  to  0.0000  N 


Figure  6-10:  Compatibility  2:  versus  Curvature  versus  Torsion 


AXIS 

X  (Curvature) 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  tenp. 
rise  of  100  deg.  C) 


SCALE 

1  unit  -  Cur  X  10.0  (Log) 
1  unit  -  Tor  X  10.0  (Log) 
1  unit  «  0.0200  B 


RAN6E 

l.OE-4  /B  to  l.OE-1  /B 
0.  1.0E“4  /B  to  l.OE-1  /B 
0.0001  B  to  0.0312  B 


U 


Compatibility  2:  U  vs  Cup  vs  Tor 


Figure  6-11:  Compatibility  2:  V  versus  Curvature  versus  Torsion 
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AXIS 

X  (Curvature) 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  tenp. 
rise  of  100  deg.  C) 


SCALE 

1  unit  >  Cur  X  10.0  (Log) 
1  unit  •  Tor  x  10.0  (Log) 
1  unit  *  0.0024  m 


RANGE 

l.OE-4  /«  to  l.OE-1  /* 

0.  l.OE-4  /«  to  l.OE-1  /m 
0.0037  «  to  0.0059  • 


Compatibility  2:  V  vs  Cur  vs  Tor 


Figure  5-12:  Compatibility  2:  V  versus  Curvature  versus  Torsion 


AXIS 

X  (Curvature) 

Y  (Toraion) 

2  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  C] 


SCALE 

1  unit  -  Cur  X  10.0  (Log) 
1  unit  -  Tor  x  10.0  (Log) 
1  unit  -  0.0005  ■ 


RAN6E 

l.OE-4  /■  to  l.OE-1  /■ 

0,  l.OE-4  /■  to  l.OE-1  /» 
0.0000  •  to  0.0005  ■ 


Compatibility  2:  W  vs  Cur  vs  Tor 


Figure  6-18:  Compatibility  2:  W  versus  Curvature  versus  Torsion 
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AXIS  SCALE 

X  (Curvature)  1  unit 

Y  (Torsion)  1  unit 

Z  (Dependent  Variable  l  unit 

plotted  for  a  temp, 
rise  of  100  deg.  C) 


RANGE 

Cur  X  10.0  (Log)  l.OE-4  /■  to  l.OE-1  /m 

Tor  X  10.0  (Log)  0.  l.OE-4  /«  to  l.OE-1  /» 

0.0200  deg  -0.0305  deg  to  0.0000  deg 

THETAX 


Compatibility  2:  Thetax  vs  Cur  vs  Tor 


Figure  6-14:  Compatibility  2: 9}^  versus  Curvature  versus  Torsion 


54 


AXIS  SCALE 

X  (Curvtture)  1  unit 

Y  (Torsion)  1  unit 

Z  (Dependent  Varleble  1  unit 

plotted  for  a  tenp. 
rise  of  100  deg.  C) 


RAN6E 

Cur  X  10.0  (Log)  l.OE-4  /•  to  l.OE-1  /■ 

Tor  X  10.0  (Log)  0.  l.OE-4  /■  to  l.OE-1  /■ 

0.0100  deg  -0.0075  deg  to  0.0000  deg 

THETAY 


Figure  6-16;  Compatibility  2:  By  versus  Curvature  versus  Torsion 
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AXIS  SCALE 

X  (CupvBture)  1  unit  -  Cup  x  10.0  (Loo) 

Y  (Topslon)  1  unit  -  Top  x  10.0  (Loq) 

Z  (Dipendont  Vspiable  l  unit  -  0.4B00  deg 

plotted  fop  a  temp. 

Plae  of  100  deg.  C) 


RANGE 

l.OE-4  /■  to  l.OE-1  /• 

0.  l.OE-4  /■  to  l.OE-1  /• 
-1.2650  deg  to  -0.0060  deg 


TOR 

Compatibility  2:  Thetaz  vs  Cur  vs  Tor 


Figure  6-16:  Compatibility  2:  versus  Curvature  versus  Torsion 


imitate  those  of  Compatibility  1;  they  show  peaks  in  the  region  where  k  =  X. 

The  forces  and  arise  out  of  the  resistance  of  the  thermal  shroud  to  the  U 
and  W  displacements  of  the  gun  tube.  Thus,  their  plots  imitate  those  of  U  and  W,  with  a 
reversal  of  sign.  Note  that  F^  is  much  smaller  than  Fj^^  in  magnitude.  This  is  because  the 
out-of-plane  IV  displacement  component  is  correspondingly  smaller  in  magnitude  than  the 
in-plane  [/  displacement,  and  as  a  result,  the  resistance  offered  by  the  shroud  to  If'  is  also 
lesser.  The  [/  and  fV  displacements  decrease  significantly  in  nnagnitude  in  comparison  with 
those  of  Compatibility  1,  as  expected,  due  to  the  restraining  forces  Fx  and  F^  resi>ectively. 
Thus,  we  find  that  in  Compatibility  2,  V  is  the  dominating  displacement  component.  The 
behavior  of  the  ffx  rotation  represents  the  straightening  of  the  helical  arc  of  Figure  5.1  and 
the  consequent  reduction  of  its  deviation  from  the  plane.  Thus  0x  remains  ^dways 

negative. 

6.4.2  Variations  on  Compatibility  2 

We  now  study  the  effect  of  variations  of  the  axial  and  flexural  rigidity  of  the 
shroud  on  the  forces  and  displacements  that  the  gun  tube  experiences.  If  either  of  the  axial 
(EA)  or  transverse  (EJ)  stiffnesses  of  the  thermal  shroud  is  increased,  as  expected,  the  force 
acting  on  the  gun  tube  increases.  If  the  shroud  is  made  more  axially  stiff,  that  is,  if  EA  is 
increased,  the  Fy  component  of  force  increases  more  significantly  than  the  Fx  and  F^ 
components.  The  reverse  is  true  of  the  case  in  which  the  shroud  is  stiffened  laterally  by 
increasing  El.  This  is  evident  from  figures  5.17, 5.18,  5.20,  5.21,  5.23  and  5.24. 

Comparing  figures  5.9  and  5.20,  we  see  that  the  increase  in  Fy  due  to  lateral 
stiffening  is  more  significant  for  the  larger  values  of  curvature.  For  mrutll  values  of 
curvature,  this  increase  is  insignificant.  This  is  because  when  the  curvature  is  small,  the 
displacement  of  the  gun  tube  end  is  primarily  along  the  direction  so  that  the  increased 
El  has  no  bearing  on  the  resistance  offered  by  the  shroud  to  the  motion  of  the  gun  tube 
end.  Thus  Fy  is  unaffected  by  increased  El  for  small  curvature.  On  the  other  hand,  for 
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Y  (Torsion) 

1  (Dependent  Variable 
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rise  of  100  deg.  C] 
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RANGE 
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0.  l.OE-4  /m  to  l.OE-1  /« 
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FORCEX 


CUfl 


TOR 


Compatibility  2. 


El  high:  Fx  vs  Cur  vs  Tor 


Figure  6-17:  Compatibility  2,  F/high:  versus  Curvature  versus  Torsion 
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AXIS  SCALE 
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Figure  5-18:  Compatibility  2,  EA  high;  Fx  versus  Curvature  versus  Torsion 
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Z  (Dependent  Variable 
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Compatibility  2.  EA  low:  Fx  vs  Cur  vs  Tor 


Figure  6-19:  Compatibility  2,  EA  low:  Fy  versus  Curvature  versus  Torsion 
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Compatibility  2.  El  high:  Fy  vs  Cur  vs  Tor 


Figure  6-20:  Compatibility  2,  £7  high;  Fy  versus  Curvature  versus  Torsion 
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Compatibility  2.  EA  high:  Fy  vs  Cur  vs  Tor 


Figure  6-21;  Compatibility  2,  EA  high:  Fy  versus  Curvature  versus  Torsion 
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AXIS  SCALE 

X  (Curvature]  1  unit 
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Figure  5-22:  Compatibility  2,  EA  low:  Fy  versus  Curvature  versus  Torsion 
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Z  (Dependent  Variable 
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Compatibility  2.  El  high:  Fz  vs  Cur  vs  Tor 


Figure  5-28;  Compatibility  2,  £/high:  versus  Curvature  versus  Torsion 


X  (Curvature)  1  unit  -  Cur  x  10.0  (Log)  l.OE-4  /■  to  l.OE-1  /■ 

Y  (Tore ion)  1  unit  -  Tor  x  10.0  (Log)  0.  l.OE-4  /n  to  l.OE-1  /■ 

Z  (Dependent  Variable  1  unit  -  12.0000  N  -1S.035B  N  to  0.0000  N 

plotted  for  a  temp, 
rise  of  100  deg.  C) 

FORCEZ 


Compatibility  2.  EA  high:  Fz  vs  Cur  vs  Tor 


Figure  6-24:  Compatibility  2,  EA  liigh:  Fr  versus  Curvature  versus  Torsion 
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Z  (Dependent  Variable 
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Compatibility  2.  EA  low:  Fz  vs  Cur  vs  Tor 


Figure  6-26:  Compatibility  2,  EA  low:  versus  Curvature  versus  Torsion 
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large  curvatures,  the  gun  tube  end  experiences  lateral  motion  to  which  the  shroud,  due  to 
its  increased  EU  offers  increased  resistance.  Thus  the  overall  resistance  increases  and  so 
does  Fy  Referring  to  Figure  5.21,  we  see  that  the  increase  in  Fy  due  to  axi£d  stiffening  is 
felt  most  for  smedl  values  of  curvature.  A  converse  reasoning  may  be  used  to  explain  this 
situation. 

If  the  EA  of  the  shroud  is  decreased,  which  amounts  to  making  it  more  easily 
extensible,  Fy  becomes  almost  constant  at  approximately  135.0  N.  This  is  explained  in 
terms  of  equations  (4.23)  as  follows:  As  the  EA  of  the  shroud  is  decreased,  we  expect  that 
the  force  Fy  decreases.  This  means  that  the  displacements  2dso  decrease,  resulting  in  lower 
values  of  Fx  and  For  the  purpose  of  this  argument,  therefore,  we  replace  the 

compatibility  equations  (4.23)  with  Equation  (4.21).  Now,  if  EA  is  decreased,  it  is  clew 
from  Equation  (4.21)  that  the  L,/(E,A^)  term  dominates  the  denominator,  so  that  Fy  is 
approximately  constant  and  equal  to  E^/ki^  T,  which  for  the  assumed  data,  evaluates  to 
147.1  N.  The  difference,  of  course,  is  due  to  the  fact  that  we  cannot  reduce  equations  (4.23) 
to  Equation  (4.21),  since  F^  and  F^  are  non-zero  in  reality.  Note  from  figures  5.19  and  5.25 
that  Fx  and  F^  are  indeed  greatly  reduced  in  magnitude  as  supposed  in  the  above 
argument. 

Another  point  worth  noting  is  that  the  graphs  of  Fx  and  do  not  show  the 
maxima  along  constant  torsion  contours  as  for  previous  compatibility  conditions.  Recall 
that  in  Compatibility  1,  we  explained  the  maxima  in  the  U  and  plots  as  representing  the 
meeting  point  of  two  conflicting  tendencies,  the  increase  in  flexibility  of  the  rod  with 
increasing  curvature,  and  the  decrease  in  force  resulting  from  the  solution  of  Equation 
(4.21).  Here,  however,  the  force  Fy  is  almost  constant,  so  that  these  maxima  are  obviously 
out  of  the  question. 

If,  on  the  other  hand,  we  decrease  the  EJ  of  the  shroud.  Compatibility  2  reduces, 
in  the  limit,  to  Compatibility  i,  since  the  shroud  becomes  more  flexible  and  offers  negligible 
resistance  to  the  lateral  motion  of  the  end  of  the  gun  tube.  The  plots  for  low  El  have  not 
been  included  for  precisely  this  reason. 
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Compatibility  2. 
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Figure  6-26:  Compatibility  2,  El  high:  U  versus  Curvature  versus  Torsion 


X  (Curvature)  1  unit  -  Cur  x  10.0  (Log)  l.OE-4  /■  to  l.OE-1  /• 

Y  (Torsion)  1  unit  •  Tor  x  10.0  (Log)  0.  l.OE-4  /■  to  l.OE-1  /■ 

2  (Dependent  Variable  1  unit  -  0.0200  m  0.0003  m  to  0.0409  m 

plotted  for  a  temp, 
rise  of  100  deg.  C) 
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Compatibility  2,  EA  high:  U  vs  Cur  vs  Tor 


Figure  6-27:  Compatibility  2,  EA  higli:  V  versus  Curvature  versus  Torsion 
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Compatidility  2. 


EA  low:  U  vs  Cur  vs  Tor 


Figure  6-28;  Compatibility  2,  EA  low:  U  versus  Curvature  versus  Torsion 
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Compatibility  2. 


El  high:  V  vs  Cur  vs  Tor 


Figure  5-29:  Compatibility  2,  FI  high:  V  versus  Curvature  versus  Torsion 


71 


AXIS 

X  (Curvature) 

Y  (Torsion) 

1  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  Cl 


SCALE 

1  unit  <■  Cur  X  10.0  (Log) 
1  unit  -  Tor  x  10.0  (Log) 
1  unit  -  0.0024  n 


RANGE 

l.OE-4  /■  to  l.OE-1  /■ 

0.  l.OE-4  /n  to  l.OE-1  /• 
0.0063  b  to  0.0063  a 


TOR 


CUR 


Compatibility  2.  EA  high:  V  vs  Cur  vs  Tor 


Figure  6-30:  Compatibility  2,  EA  high;  V  versus  Curvature  versus  Torsion 
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Compatibility  2. 


EA  low;  V  vs  Cur  vs  Tor 


Figure  6-31:  Compatibility  2,  EA  low:  V  versus  Curvature  versus  Torsion 
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Compatibility  2.  El  high:  W  vs  Cur  vs  Tor 


Figure  B-82:  Compatibility  2,  ^7  high:  W  versus  Curvature  versus  Torsion 
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Compatibility  2.  EA  high:  W  vs  Cur  vs  Tor 


Figure  5>83:  Compatibility  2,  EA  high:  W  versus  Curvature  versus  Torsion 
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Figure  6>84:  Compatibility  2,  EA  low:  W  versus  Curvature  versus  Torsion 
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We  observe  from  figure  5.26  and  5.32  that  when  the  El  of  the  shroud  is  increased, 
the  transverse  displticements,  U  and  W,  decrease  considerably,  although  there  is  no  change 
in  the  shape  of  the  plots  as  compared  to  those  of  Compatibility  2.  The  V  displacement 
shown  in  Figure  5.29  decreases,  for  increased  El,  noticeably  for  higher  values  of  curvature, 
from  0.0059  m  in  Compatibility  2  (Figure  5.12),  to  0.0057  m  in  the  present  case,  since  Fy 
increases  for  higher  values  of  curvature  as  explained  earlier. 

If  EA  of  the  shroud  is  increased,  the  V  displacement  attains  an  almost  constant 
value  of  0.0063  m.  This  is  shown  in  Figure  5.30.  This  is  because  the  shroud  has  now 
become  almost  infinitely  rigid  as  comp£Lred  to  the  gun  tube.  Thus  the  curvature  and 
torsion  of  the  gun  tube  are  irrelevant.  The  magnitude  of  V  is,  in  the  limit,  simply  equal  to 
the  expansion  of  the  thermal  shroud,  that  is,  X/tA  T.  The  transverse  displacements  U  and 
W  increase  from  their  values  in  Compatibility  2.  This  may  be  attributed  to  the  enormous 
increase  in  the  force  exerted  by  the  shroud  on  the  gun  tube. 

If  the  EA  of  the  shroud  is  decreased,  we  find  that  the  maxima  along  lines  of 
constant  torsion  seen  in  the  plots  of  U  and  W,  and  the  inflexion  point  in  the  plots  of  V 
disappear.  This  is  due  to  the  fact  that  the  force  exerted  by  the  shroud  on  the  gun  tube  is 
almost  a  constant,  as  can  be  seen  from  Figure  5.22.  All  of  the  above  displacements  decrease 
considerably  in  magnitude,  mainly  due  to  the  fact  that  the  force  exerted  by  the  shroud  on 
the  gun  tube  is  very  small. 

The  rotations  9x'  ^  Yi  are  the  angular  displacements  of  the  tangent  to  the 

centroideil  axis  of  the  gun  tube.  These  therefore  represent  the  straightening  of  the  helical 
arc  as  it  is  extended.  The  plots  of  6^,  ^y,  and  Og  of  figures  5.35  -  5.43  are  generally  similar 
to  those  of  Compatibility  2.  If  either  of  the  El  or  EA  of  the  shroud  is  increased,  all  of  the 
rotations  increase  since  the  gun  tube  tends  to  straighten  to  a  greater  extent,  on  account  of 
the  shroud’s  increased  stif&ess.  If  the  EA  of  the  shroud  decreases,  then  the  maxima  of  the 
plots  of  Ox  and  0^  that  occur  in  the  region  x  =  A,  disappear.  This  is  because  the  force 
exerted  on  the  gun  tube  by  the  shroud  approaches  a  constant  value. 
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Compatibility  2.  El  high;  Thetax  vs  Cur  vs  Tor 


Figure  5-SB:  Compatibility  2,  El  high:  Ov  versus  Curvature  versus  Torsion 
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Compatibility  2.  EA  high:  Thetax  vs  Cur  vs  Tor 


Figure  6*36:  Compatibility  2,  EA  high:  versus  Curvature  versus  Torsion 
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Compatibility  2,  EA  low:  Thetax  vs  Cur  vs  Tor 


Figure  6-87:  Compatibility  2.  EA  low:  versus  Curvature  versus  Torsion 
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Figure  6-88:  Compatibility  2,  F/high:  8 y  versus  Curvature  versus  Torsion 
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Compatibility  2.  EA  high:  Thetay  vs  Cur  vs  Tor 


Figure  6-39:  Compatibility  2,  EA  high:  0yr  versus  Curvature  versus  Torsion 


SCALE  RANGE 

1  unit  -  Cur  X  10.0  (Log)  l.OE-4  /•  to  l.OE-1  /• 

1  unit  -  Tor  x  10.0  (Log)  0.  l.OE-4  /■  to  l.OE-1  /• 

1  unit  -  0.0100  deg  -O.OOBO  deg  to  0.0000  deg 


THETAY 


X  (Curvature)  1  unit  *  Cur  x  10.0  (Log) 

Y  (Toraion)  1  unit  >  Tor  x  10.0  (Log) 

Z  (Dependent  Variable  i  unit  «  0.0001  deg 

plotted  for  a  temp. 

rlae  of  100  deg.  C)  THETAY 


Compatibility  2.  EA  low:  Thetay  vs  Cur  vs  Tor 


l.OE-4  /■  to  l.OE-1  /■ 

0,  l.OE-4  /m  to  l.OE-1  /■ 
-O.OOOi  deg  to  0.0000  deg 


Figure  5-40:  Compatibility  2,  EA  low:  9  y  versus  Curvature  versus  Torsion 
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AXIS  SCALE  RANGE 

X  (Curvature)  1  unit  -  Cur  x  10.0  (Log)  l.OE-4  /■  to  l.OE-1  /« 

Y  (Torsion)  1  unit  -  Tor  x  10.0  (Log)  0,  l.OE-4  /«  to  l.OE-1  /• 

Z  (Dependent  Variable  i  unit  -  0.4B00  deg  -0.BB46  deg  to  -0.0036  deg 

plotted  for  a  temp, 
rise  of  100  deg.  C) 


TOR 

Compatibility  2.  El  high:  Thetaz  vs  Cur  vs  Tor 


Figure  6-41:  Compatibility  2,  F/high:  9g  versus  Curvature  versus  Torsion 
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AXIS 

X  (Curvature) 

Y  (Torsion] 

Z  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  C) 


SCALE 

1  unit  -  Cur  X  10.0  (Log) 
1  unit  -  Tor  x  10.0  (Log) 
1  unit  -  0.4B00  deg 


RANGE 

l.OE-4  /■  to  l.OE-1  /m 
0.  l.OE-4  /*  to  l.OE-1  /* 
>1.6577  deg  to  -0.0102  deg 


CUR 


Compatibility  2. 


TOR 

EA  high:  Thetaz  vs  Cur  vs  Tor 


Figure  5*42:  Compatibility  2,  EA  high:  versus  Curvature  versus  Torsion 


X  (Curvature) 

Y  (Torsion) 

Z  (Dependent  Variable 
plotted  for  a  temp, 
rise  of  100  deg.  C) 


i  unit  •  Cur  X  10.0  (Log) 
1  unit  >  Tor  X  10.0  (Log) 
1  unit  -  0.0100  deg 

THETAZ 


l.OE-4  /•  to  l.OE-1  /■ 

0,  l.OE-4  /•  to  l.OE-1  /• 
-0.0140  deg  to  0.0000  deg 


CUR 


TOR 


Compability  2.  EA  low:  Thetaz  vs  Cur  vs  Tor 


Figure  6>4S:  Compatibility  2,  EA  low;  07  versus  Curvature  versus  Torsion 
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The  following  chapter  contains  the  formulation  of  a  general  space-curved  finite 
element  which  is  later  used  to  independently  confirm  the  values  of  the  more  exact 
displacements  found  from  the  first  method  of  solution. 


CHAPTER  6 


A  GENERAL  SPACE-CURVED  FINITE  BEAM  ELEMENT 


6.1  Introduction 

In  the  previous  chapters,  we  formulated  and  performed  a  parametric  study  of  an 
exact  solution  to  the  problem  of  finding  the  displacements  of  the  end  of  the  gun  tube.  In 
the  present  chapter,  we  present  an  alternative  approach:  the  development  of  a  general 
space-curved  finite  beam  element  within  which  variable  curvature  and  torsion  are 
I  permitted.  By  the  use  of  an  assemblage  of  such  elements,  we  may  model  a  rod  of  arbitrarily 

varying  cross-section,  and  curvature  and  torsion. 

The  element  is  derived  by  employing  the  variational  method.  The  formulation 
assumes  the  expression  for  the  deformation  str^  energy  density  of  a  curved  beam,  whose 
centroidal  axis  is  an  arbitrary  space  curve  of  variable  torsion  and  curvature,  obtained  by 
Tsay  and  Kingsbury.  [Tsay  86]  This  expression  is  derived  from  the  linearized  strain- 
displacement  relations  due  to  Kingsbury  [Kingsbury  84],  which  assume  that  cross-sections 
remain  plane  and  normal  to  the  centroidal  curve  under  deformation. 

The  most  important  feature  of  this  element  is  its  generality.  The  element  sti&ess 
matrix  obtained  is  valid  for  any  beam  element  whose  centroidal  axis  is  a  curve  with  a 
constant  rate  of  change  of  curvature  and  tornon  throughout  its  length  and  also  a  constant 
shape  of  cross-section.  In  the  following  chapter,  we  shall  present  a  comparison  between  the 
exact  and  the  finite  element  solutions. 

( 
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6.2  The  Displacement  Field 

6.2.1  Assumptions  of  the  Formulation 

The  differential  equations  describing  the  motion  of  space-curved  beams  have  been 
formulated  by  Kingsbury  in  terms  of  the  displacement  components  v,  v,  w  and  4>  which 
represent  four  of  the  six  unknown  displacements  or  degrees  of  freedom  permitted  to  any 
point  on  the  beam.  Exact  solution  of  these  four  differential  equations  is  made  difficult  by 
the  fact  that  they  are  involve  coupling  in  the  above  four  variables  due  to  the  introduction  of 
torsion,  u,  v  and  w  are  translational  displacements  of  the  centroid  of  a  generic  cross-section 
in  the  directions  of  the  local  principal  normal,  binormal  and  tangent  respectively,  and  (j> 
represents  the  angular  displacement  of  any  point  of  the  cross-section  about  the  local 
tangent.  These  four  displacement  quantities  vary  along  the  length  of  the  beam;  this  may  be 
expressed  as 

li  =  u{z) 

V  =  v{r) 
w  =  w{z) 

<t>  =  0{z)  (6.1) 

where  ^  is  the  arc  length  measured  along  the  space-curved  beam.  The  remaining  two 
angular  displacements  sire  derived  fKingsbury  84]  from  u(z),  v{z)  and  w(r),  according  to  the 
equations  below,  which  follow  from  the  neglect  of  transverse  shear  deformation: 


dv  . 

Q,  =  -  —  -An 
*  dz 

(6.2) 

du  . 

0„  =  -  -  Xv  +  KV! 

*  dz 

(6.3) 

The  second  assumption  made  is  that  angular  displacements  are  small  and  hence  may  be 
transformed  from  one  coordinate  system  to  another  the  same  way  as  translational 
displacements  are. 
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6.2.2  The  Displacement  Field:  A  Polynomial  Representation 


Consider  the  finite  element  shown  in  Figure  6.1.  Each  node  of  the  element  has  six 
degrees  of  freedom,  three  translational  emd  three  rotational.  The  element,  therefore,  has 
twelve  degrees  of  freedom  in  all.  As  stated  in  the  foregoing  section,  the  four  displacements 
in  terms  of  which  the  differential  equations  are  formulated  ue  functions  of  the  independent 
variable  z.  the  arc  length  along  the  beam. 

Noting  that  the  highest  orders  of  derivatives  of  these  displacements,  u,  v.  vi  and 
appearing  in  the  four  coupled  differential  equations  of  Kingsbury  are  four,  four,  two  and 
two  respectively,  we  adopt  the  following  polynomial  representation  of  the  displacement 
field: 


u{z)  =  —  a2Z—  —  a^z' 


v{z)  =  —  b^z—  b^z^  ~  b^z^ 


iy(r)  =  cj  —  Cjz 


o{z)  =  —  djz  (6.4) 

The  reader  may  note  that  the  total  number  of  undetermined  constants  in  the  above 
equations  is  the  same  as  the  number  of  element  degrees  of  freedom  to  be  specified,  that  is, 
twelve.  It  may  be  mentioned  here  that  the  above  displacement  field  satisfies  the  differential 
equations  of  Kingsbury  exactly  in  the  case  that  curvature  k  and  torsion  A  (and  their 
derivatives)  are  identically  zero,  in  other  words,  in  the  straight  beam  case.  It  does  not, 
however,  satisfy  the  differential  equations  for  the  space-curved  beam.  However,  in 
consideration  of  the  fact  that  we  are  dealing  with  small  curvatures  and  torsions 
representing  the  imperfections  in  the  straightness  of  gun  tubes,  this  assumption  appears  to 
be  reasonable. 
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x(L),  u2 


t(L).  w2 


x(0),  ul 


z(0).  wl 


y(0),  vl 


y(L).  v2 


X  s  Principal  Nonna! 
y  s  Binonnal 
z  «  Tangent 


Figure  6-1:  Degrees  of  Freedom  of  the  Finite  Element 

Having  chosen  polynomials  to  describe  the  displacements,  we  now  seek  to  obtain 
the  shape  functions  for  the  displacement  field.  To  achieve  this,  we  must  first  determine  the 
unknown  constants  introduced  in  equations  (6.4).  These  constants  are  expressed  in  terms 
of  the  values  of  the  element  d.o.f.  (degrees  of  freedom),  which  are  specified.  Differentiating 
equations  (6.4)a  and  (6.4)b,  we  have: 


U  (z)  =  —  =  Oj  i-  203Z-r  3a4Z' 

dz 


v'{z)  =.^  =  6}  4-  263Z+  364Z^ 
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Let  Uj,  V|,  Uj <j)i,  and  Uj,  Vj,  tU})  U} ',  Vj ^2  values  of  the  element  d.o.f.  at 

the  two  nodes  given  by  z  =  0  and  z  =  I  respectively.  Eveduating  equations  (6.4)  and  (6.5)  at 
the  two  nodes,  setting  them  equal  to  the  known  values  of  the  d.o.f.  and  inverting  the 
relationships,  we  arrive  at  the  following  expressions  for  the  constants: 

fli  =  txi 


“2-^1 


_  3(tt2-  Ui)  _  (2lii'  -r  ttj') 


/ 

2(«,- 

“2)  , 

_  (“1 

'  -h') 

“4  - 

61  = 

n 

62  = 

63  = 

3  (u2  “ 

~T' 

111  - 

.  (2vi 

■  -  V2  ') 

1 

64  = 

2(vi- 

J^2)  ^ 

-Vj') 

Cl  =  Wi 


,  _  (Wj  -  Wi) 

’  r~ 

di  =  <?! 

-  {4>2~  ^l) 
’  I 


(6.6) 
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6.2.3  The  Shape  Function  Representation 


Having  determined  the  unknown  constants  in  terms  of  the  known  nodal  values  of 
the  displacements,  we  may  now  write  the  shape  functions  for  the  assumed  displacement 
field.  Substituting  equations  (6.6)  into  equations  (6.4)  and  rearranging  the  terms,  we 
obtain  the  shape  functions  for  the  displacements  u  and  v  to  be  the  following  cubic 
polynomials  in  z. 


I,(r)  =  l- 


L2{z) 


z 


z^ 


hiz)  = 


2z^ 


Uz)  =  - 


We  also  obtmn  the  following  linear  shape  functions  for  displacements  w  and 


(6.7) 


iVi(r)  =  l-5 


(6.8) 


Combining  the  shape  functions  of  equations  (6.7)  and  (6.8)  with  the  nodal  values 
of  the  displacements,  we  obtain  the  following  expressions  for  the  four  displacements; 


(\  M  3z*  ,  2r*  r*. 


1*  T -T’ 


+  ( 


3z*  2r’  -* 
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/  \  2z\  .  _2z*  z’. 


1 


-  ( 


Sz*  2z^ 


I*  I 


z  z* 

)  ’'2  +  ( “  “7  75  )  ’'j 


/  /- 


wiz)  =  (1-pw,  ^  (-^)«;2 


=  (1-^)  <^1  -  (p  <^2 


(6.9) 


Having  obtained  the  displacements  as  interpolations  of  their  nodal  values  in  terms 
of  the  arc  length  z,  we  are  now  in  a  position  to  get  an  expression  for  the  strain  energy  of 
deformation  of  the  beam  element  in  terms  of  the  nodal  values  of  the  displacements.  This  is 
accomplished  in  the  following  section. 


6.3  Strain  Energy  of  the  Beam  Element 

6.3.1  The  Linear  Strain  Energy  Density  Function 


The  expression  for  the  linear  strain  energy  density  (strain  energy  per  unit  arc 
length)  of  the  deformed  beam,  as  obtained  by  Tsay  and  Kingsbury  (Tsay  86j,  is: 

F  =  /_  {  F  :  - - — r  .i - —  - - — !-wv  A  -^uv  A 

2  2  2 

-t-tt«'AA'-v''K<i)-uA'K<^-u'Ax^]  T  G\  -f  ti  v'A  k  * 

2  2  2  2 

-r  7^  {  E[  v"v)k'  +  v'  'vi'k  -  V  v"A  '  -  i;'»"A  +  «"«" 

+  «ti>A'x'  +  aw'A'/c-t»»(A')*-«»'AA'  +  tttt"A'  +  ii'ti)A»c' 

+  tt'ti;'Ax  -  u'vAA'  -  tt'»'A*  +  tt'tt"A-  v)kk'<^-  w'xV 


d4 

^vX'K(t>  —  v'XK.<t>  -  u"»c<^  j  +  G[  v'wk^-tuv/Xk^  +  wk^^' 

-  vv'Xk^-*-u'v'k^~v'Xk<I>  -  11t>  (A«)  *  -  vX  K  uu'X 

—  u'kO—  UK  x^<t)  j} 

yy  2  2  2  2  2 


-  u'  v'A  -  u''vX'-!-u''wk'^u"w'k  +  vv'XX'-  v'wXk' 

—  v’w'X  K  —  vwX'k'  -  vw'X'k  +  ww'k  k'  1  +  G 


VJ^K*  3 

- - vwXk 


3  \  2^  (vXk)^  ,  .  2  \2  .  (“  'c) 

~  U  WK  —  W  X  K  - - -  -  U  vX  K  -  vX  K(p  — 

-uXK<,-i^-i^  ) 

2  2 

-A  {E  - — !—-uwk~- - 

2  2 


(6.10) 


6.8.2  Strain  Energy  of  the  Element 

Substituting  equations  (6.9)  into  Equation  (6.10),  the  strain  energy  density 
function,  and  integrating  the  result  over  the  length  of  the  element  /,  we  arrive  at  the 
expression  for  its  strain  energy  of  deformation.  This  expression,  however,  is  presented  in 
Appendix  B.  owing  to  its  length. 


0.4  The  Element  Stifhiess  Matrix 


Next,  we  apply  Castigliano’s  First  Theorem  to  the  strain  energy  expression.  The 
theorem  states  that  if  a  linearly  elastic  structure  is  subjected  to  a  displacement,  the  load  in 
its  direction  that  causes  the  displacement  is  equal  to  the  partial  derivative  of  the  strain 
energy  of  deformation  of  the  structure  with  respect  to  that  displacement.  We,  therefore, 
differentiate  the  strain  energy  with  respect  to  the  nodal  values  of  displacements  which 
appear  in  the  expression  and  thus  obtain  the  components  of  the  internal  force  and  moment 
acting  on  the  ends  of  the  element  in  their  respective  local  coordinate  systems  specified  by 
the  tangent,  binormal  and  principal  normal  at  those  points. 


95 


The  coefficients  of  the  nodal  values  of  displacements  appearing  in  the  expressions 
for  nodal  forces  and  moments  obtained  by  differentiating  the  strain  energy  are  in  fact  the 
elements  of  the  element  stiffness  matrix  K.  The  relationship  between  the  nodzd  forces  and 
the  noHai  displacements  in  terms  of  the  element  stif&iess  matrix  K  is: 

F  =  KZl.  (6.11) 

where,  F  is  the  internal  force  vector  given  by 

F  =  Vzi  Vj/i  Vz^  A/zi  Mvi  Mz^  Vx^  Vy^  Vz^  Jl/zj  H/y,  ,  (6.12) 

and  is  the  nodtil  displacement  vector  given  by 

=  til  Vj  Wi  <^xl  “vl  ^2  ’'2  ^"2  “i2  “y2  .  (6-13) 

The  a  and  a  values  in  the  nodal  displacement  vector  are  as  given  by  equations 

X  y 

(6.2)  and  (6.3).  Since  we  are  interested  in  displacements  occurring  for  small  values  of 
curvature  k  emd  torsion  A,  we  drop  the  last  term  from  Equation  (6.2)  and  the  last  two 
terms  from  Equation  (6.3).  thereby  getting  the  simpler  results 

Q*  =  -  and  a  =  ^  (6.14) 

az  *  dz 

The  above  relationships  hold  exactly  for  the  case  of  the  straight  beam,  and  are,  therefore, 
consistent  with  the  shape  functions  assumed  for  the  element. 


Note  that  this  simplification  results  in  {in  uncoupling  of  displacements,  thereby 
rendering  the  problem  more  amenable  to  solution.  We  now  have  a,i  =  -  vj ',  =  «i  \ 

~  V]'  and  =  Uj',  as  per  Elquation  (6.14),  so  that  the  moments  Mxy,  Mx^ 

and  A/V}  are  obtained  by  applying  Castigliano's  First  Theorem  as  follows: 


Afzi 


dU 

dvi ' 
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My,  = 


au 

du, ' 


Mx2  =  ~ 


au 

avj ' 


Myj  =  ^  (6.15) 

OU2 

The  term  U  appearing  in  the  above  equations  is  the  strain  energy  of  deformation  of  the 
element. 


The  12  X  12  element  stiffness  matrix  is  symmetric,  as  expected.  The  expressions 
for  the  elements  of  the  stiffness  matrix  in  terms  of  the  geometrical  parameters  defining  the 
element,  k,  k  '.  A,  A  ',  /,  cross-sectional  properties  lyy,  A,  and  material  properties  E 

and  are  provided  in  Appendix  D,  which  contains  a  finite  element  program  used  to 
compeu’e  the  displacements  obtained  by  using  the  general  space-curved  beam  element  with 
those  of  the  exact  solution  presented  in  previous  chapters. 

6.4.1  Reduction  to  the  Straight  Beam  Case 

In  this  subsection,  we  present  a  check  on  the  element  stifihess  matrix  obtained 
earlier  for  the  case  of  the  straight  beam  element.  For  the  straight  beam,  we  have  zero 
torsion  and  curvature.  Substituting  these  values  into  the  expressions  found  for  the 
elements  of  the  stiffness  matrix,  we  obtain  the  following  results  which  agree  with  the 
standard  straight  beam  element  matrix  [Zienkiewicz  77;.  (Only  the  non-zero  upper 
triangular  elements  are  listed  row-wise,  since  the  matrix  is  symmetric.) 

=  -  K{1,7) 

K[l,2)  =  ~ 

K{IA)  =  if(l,10)  =  - 


K(1,5)  =  JC(1.11)=^^ 


K{2,2)  =  -  K{2.S)  = 


K(2A)  =  A-(2,10)  =  - 


/f(2,5)  =  K(2,11) 
^•(2,7)  =  - 


^(3.3)  =  -  ^-(3,9)  = 


£A 


A,  El 

K{AA)  =  2if(4.10) 


if(4,5)  =  2K{AX\)  -  ~ 
if(4J)  =  -i^ 
if(4,8)  =  - 


/f(5,5)  =  2/f(5,ll) 
/C(5,7)  =  - 
if(5,8)  =  - 


K{hAO)  =  - 


/r(6,6)  =  -  Jf(6,12)  =  — 
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K{7J)  = 
K{7.S)  = 
/f(7,10)  = 
/C(7,ll)  = 
K{%.8)  = 
/f(8.10)  = 
/f(8,ll)  = 
/f(9,9)  = 
/C(10,10) 

1^(12,12) 


12  El. 


mi 


12  El. 


SI 


QEL 


SI 


6  El 


sa 


12^4 


QEI.^ 


&El 


SI 


/ 


iEL 


4Bl 


SI 


4El 


mi 


(6.16) 
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6.6  A  Finite  Element  Example  Solution 

In  this  section,  we  present  the  results  obtained  by  using  the  finite  element 
program  of  .4ppendi.'<  D  for  the  displacements  of  the  end  of  the  gun  tube,  whose  model  is 
shown  in  Figure  6.2.  The  gun  tube  is  modeled  as  a  cantilevered  beam  with  three  elements. 
Note  that  compatibility  of  displacements  of  the  ends  of  the  gun  tube  and  the  shroud  is  not 
enforced  but  we  use  as  an  approximation,  the  value  of  the  force  found  for  the  exact 
solution,  for  the  corresponding  values  of  initial  curvature  and  torsion.  The  shroud  is 
removed,  but  the  approximate  shroud  force  is  applied  to  the  gun  tube  in  a  direction  parallel 
to  the  chord  joining  the  ends  of  the  gun  tube. 


Fixed  End 


Free  End 


14  cm 


1 

jecm 


10  cm 


1.75  m 
Element  1 


1.75  m 
Element  2 


1.75  a 
Element? 


Figure  6-2:  Model  of  a  Gun  Tube 

It  must  be  brought  to  attention  that  the  global  stiffness  equation  whose  solution  is 
implemented  by  the  finite  element  program  is  formulated  in  the  thori.  coordinate  system. 
Since  the  element  stiffness  matrix  is  developed  in  the  local  tangent  system,  we  tnake  use  of 
the  transformations  of  Section  4.3  in  obtaining  the  sti&ess  matrix  in  the  chord  system. 
This,  although  not  a  necessity,  facilitates  the  comparison  between  the  displacements  found 
from  the  exact  and  the  finite  element  solutions  which  is  presented  in  the  following  chapter. 


( 
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The  solution  of  the  global  stiSness  equation  in  the  chord  system,  F  =  KA, 
proceeds  as  follows:  We  first  partition  the  global  force  and  displacement  vectors  into  blocks 
of  known  (denoted  by  the  subscript  a)  and  unknown  (denoted  by  the  subscript  0) 
quantities,  and  rearrange  the  stifhiess  matrix  accordingly.  Thus  the  global  equation  may  be 
written  as: 


(0  ■  (&fc)  (0 

Now,  since  one  end  of  the  gun  tube  is  held  fixed,  the  known  displacements  are  zero. 
Thus  the  unknown  displacements  and  forces  (reactions)  may  be  obtained  from  Equation 
(6.17)  as: 


-3  = 

Fa  =  K3,K„„-‘F„  (6.18) 


The  general  finite  element  allows  for  variation  of  the  curvature  and  torsion  along 
its  length.  In  the  present  example,  however,  this  facility  is  not  used.  Data  and  results  for 
this  example  follow. 


Global  Data: 


Number  of  elements  =  3 
Length  of  gun  tube  =  5.25  m 

Curvature  k  =  1.0  x  10”^  per  m 
Torsion  A  =  1.0  xl0~®  perm 
Rate  of  change  of  curvature  x '  =  0 
Rate  of  change  of  torsion  A '  =0 

Modulus  of  elasticity  F  =  4.35097  x  10*  N/m* 
Shear  modulus  G  =  1.67345  x  10*  N/m* 

Force  exerted  by  thermal 
shroud  along  the  chord  = 


61441.952  N 
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Elemental  Data: 


Length  of  each  element  = 

1.75  m 

Cross-section  of  element  1: 

Inner  radius 

=  0.06 

m  = 

6  cm 

Outer  radius 

=  0.14 

m  = 

14  cm 

Cross-sections  of  elements  2  and  3: 

Inner  radius 

=  0.06 

m  = 

6  cm 

Outer  radius 

=  0.10 

m  = 

10  cm 

Results: 

Computed  Reactions  at  the  Fixed  End  in  the  Chord  System: 


^Xc 

= 

7.409  X  10"®N 

= 

-61441.952  N 

= 

1.964  X  10"^N 

-X. 

= 

-7.409  X  10"^  N-m 

= 

-1.954  X  10"^  N-m 

= 

-4.439  X  10"^  N-m 

End  Displacement  of  the  Gun  Tube  in  the  Chord  System: 


AX  = 

5.984  X  I0"'‘m 

Ay  = 

2.950  X  10"’ m 

AZ  = 

1.510  X  10"®m 

A^jf  = 

3.712  X  10"’  deg 

= 

3.266  X  10"*  deg 

II 

o 

-1.146  X  10" ’deg 

Note  that  the  applied  force  acts  along  the  chord  and  so  passes  through  the  point 
at  which  the  gun  tube  is  encastred.  Thus  the  moments  generated  at  the  fixed  end  must  be 
zero;  so  should  the  reactions  in  directions  perpendicular  to  the  chord.  This  is  borne  out  by 
the  computed  reactions  at  the  fixed  end  shown  above,  of  course,  with  some  error.  In  the 
next  chapter,  we  shall  present  a  comparison  of  the  exact  and  finite  element  solutions. 
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CHAPTER  T 

A  COMPARISON  OF  THE  EXACT  AND  FINITE  ELEMENT  ANALYSES 

7.1  The  Impact  of  Simplifying  Assumptiona 

Recall  that  the  displacement  field  that  we  assumed  for  the  finite  element 
formulation  of  Chapter  6  satisfied  the  differential  equations  of  Kingsbury  [Kingsbury  84] 
only  in  the  case  of  the  straight  beam.  In  other  words,  if  we  were  to  substitute  these 
displacements  and  their  derivatives  into  the  differential  equations,  we  would  expect  a  non¬ 
zero  residual.  This  implies  a  certain  inherent  inaccuracy  in  the  finite  element,  and  we 
naturally  expect  that  the  element  will  show  a  greater  accureu:y  for  the  smaller  values  of 
curvature  and  torsion. 

It  must  also  be  noted  that  the  simplification  of  Elquation  (6.14)  would  affect  the 
accuracy  of  all  the  displacements,  particularly  w,  a,  and  q^. 

7.2  A  Comparison  of  the  Exact  and  Finite  Element  Methods 

Presented,  in  this  section,  are  some  of  the  displacements  obt^ed  by  the  exact 
method  under  the  first  compatibility  condition  described  in  Section  4.5,  and  the 
corresponding  displacements  obtained  using  the  general  beam  element  of  Chapter  6,  for  a 
few  selected  values  of  initial  curvature  and  torsion  of  the  gun  tube.  These  displacements 
were  computed  using  the  programs  given  in  appendices  C  and  D  respectively. 

The  applied  forces  used  in  the  finite  element  program  were  supplied  from  the 
output  of  the  program  for  generating  the  exact  solutions,  in  which  compatibility  conditions 
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were  enforced.  The  gun  tube  and  shroud  data  assumed  for  the  comparison  were  the  same 
as  those  given  in  Chapter  5,  and  for  ease  of  reference,  they  are  presented  below; 

Gun  Tube: 


Length  =  5.25  m 

Inner  radius  =  0.06  m 

Outer  radius  =  0.10  m 

Elasticity  modulus  =  4.35097  x  10®  N/m* 

Shear  modulus  =  1.67345  x  10  ®  N/m* 

Thermal  Shroud: 

Inner  radius  =  0.11  m 

Outer  radius  =  0.145  m 

Temperature  rise  =  100  °  C 

Elasticity  modulus  =  4.35097  x  10  ®  N/m* 

Coefficient  of  thermal  expansion  =  1.206  x  10  “  ®  /  °  C 


The  finite  element  analysis  of  the  gun  tube  was  performed  with  a  3-element  tmd  a 
10-element  model  to  verify  the  expected  increase  in  the  accuracy  of  the  solution  against  the 
exact  solution.  Five  different  cases  of  initial  curvature  and  torsion  are  presented  in  this 
chapter  with  conclusions  about  the  results  inserted  at  suitable  junctures.  Note  that  the 
exact  reactions  at  the  fixed  end  of  the  gun  tube  consist  of  only  one  equal  but  opposite  force 
-  Fyc-  Since  the  force  due  to  the  shroud  acts  along  the  chord  itself,  no  reaction  moments 
are  produced  at  the  fixed  end. 

CASE  1: 

K  =  1.0  X  10  “  ^  /  m 
A  =  0.0  /m 

Applied  force  in  the  chord  system,  Fy^  =  61441.9292446  N 

Computed  Reactions  in  the  Chord  System,  in  N  and  N-m: 

3  Elements  10  Elements 

F^^  -  5.2988725  x  10 "  7.2120088  x  10 '  ” 

Fy^  -61441.9292446  -  61441.9292446 
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^Zc 

0.0 

0.0 

^Xc 

0.0 

0.0 

^Yc 

0.0 

0.0 

^Zc 

-1.4694506x10’^ 

-1.3903836x10’* 

Displacements  in  the  Chord  System, 

in  m  and  deg: 

Exact  Solution 

3  Elements 

AX 

C 

6.5259152  X  10 

6.5387099  X  10 

Ay 

c 

3.6873859x10’* 

3.6873662x10’* 

AJ 

c 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

^9zc 

-1.4272158  x  10’* 

-1.4293170x10’* 

CASE  2: 

K  =  1.0X10’*  'm 

A  =  0.0  / 

m 

Applied  force  in  the  chord  system,  F 

=  55426.4228976  N 

Computed  Reactions  in  the  Chord  System,  in  N  and  N-m: 

3  Elements 

10  Elements 

F 

Xc 

-4.4158832x10’* 

-4.4949355x10’* 

F 

Yc 

-55426.4279317 

-55426.4294723 

F 

Ze 

0.0 

0.0 

^Xc 

0.0 

0.0 

^Yc 

0.0 

0.0 

-121.4720064 

- 12.4268252 

Displacements  in  the  Chord  System, 

in  m  and  deg: 

Exact  Solution 

3  Elements 

<I 

5.8975078x10’* 

5.4112829x10’* 

Ay 

3.9458062  X  10’ ’ 

3.9069921x10’* 

AZ 

0.0 

0.0 

C 


10  Elements 

6.5387628  x  lO^* 
3.6873663  X 10 
0.0 
0.0 
0.0 

-1.4293262x10’* 


10  Elements 

5.8475782x10’* 

3.9410363x10’* 

0.0 
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0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

- 1.2873947 

- 1.2035964 

-1.2799743 

From  the  above  figures  for  cases  1  and  2,  we  observe  that  the  finite  element 
method  yields  very  accurate  values  of  displacements,  for  zero  torsion,  when  the  curvature  is 
small.  As  mentioned  following  the  finite  element  example  of  the  previous  chapter,  an 
indication  of  the  accuracy  of  the  solution  is  the  accuracy  of  the  computed  reactions.  We 
expect  that  the  only  non-zero  reaction  at  the  fixed  end  of  the  gun  tube  is  an  opposite  force 
to  that  exerted  by  the  thermal  shroud. 


With  this  in  mind,  it  is  easy  to  see  that  in  the  second  case,  the  accuracy  of  the 
solutions  is  not  as  good  as  the  first.  However,  the  error  in  the  moment  M  ,  whose  value 
should  be  zero,  is  still  not  very  large  considering  the  enormous  magnitude  of  the  applied 
force.  Note  that  in  Case  1.  the  10-element  model  produces  slightly  more  accurate  reactions 
than  the  3-element  model.  There  is  also  not  much  improvement  in  displacement  values 
from  the  3-element  model  to  the  10-element  model.  This  lack  of  improvement  may  be 
attributed  to  the  inaccuracy  of  the  element  itself,  mberent  in  the  approximate  displacement 
field  assumed  in  its  derivation.  In  Case  2,  however,  the  improvement  is  much  more  visible. 

CASE  3: 

K  =  1.0  X  10“^  /m 
A  =  1.0x10"'*  /m 

Applied  force  in  the  chord  system,  =  61441.9460865  N 

Computed  Reactions  in  the  Chord  System,  in  N  and  N-m: 

10  Elements 

7.4090253x10"^ 

-61441.9460865 
2.6942148x10"^° 

-7.4090355x10"^ 


I 


3  Elements 

7.4089804x10"^ 
f’y*  -61441.9460865 

F,  2.6942000x10"*° 

Zc 

-7.4090486x10"’ 

AC 
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"yc 

-7.4802397x10"^ 

-7.4802703x10"^ 

- 1.4889072  X  10 

-1.5849468x10"® 

Displacements  in  the  Chord  System, 

in  m  and  deg: 

Exact  Solution 

3  Elements 

10  Elements 

AX 

c 

6.5390314  X  10 

6.5387117x10"'* 

6.5387646  X 10  "■* 

<] 

3.6873852x10"^ 

3.6873672x10"® 

3.6873674X10"® 

<1 

1.1735696x10"® 

1.4904136x10"^ 

1.4875948x10"^ 

8.6536987x10"^^ 

3.7464527x10"* 

3.7464405x10"® 

-4.6529785x10"“ 

2.3128574x10"® 

2.3147083x10"® 

-1.4272162x10"* 

-1.4293173x10"* 

-1.4293266x10"* 

It  is  important  to  observe  in  Case  3  that  even  though  the  reactions  have  been 
computed  very  2u:curately  (as  compared  to  Case  2,  for  example),  the  displacements  ^2^, 
which  are  ’’out-of- plane”  in  the  sense  of  Section  5.2  show  fairly  large  errors  for 
both  the  3-element  and  the  10-element  models.  This  is  because  the  introduction  of  torsion 
results  in  a  coupling  of  displacements  which  is  not  accounted  for  in  the  finite  element 
model.  Thus  the  finite  element  solution  is  unresponsive,  in  terms  of  its  accuracy  with 
regeird  to  these  displacements,  to  an  increase  in  the  number  of  elements. 


The  assumptions  that  we  made  in  Chapter  6  yield  simplified  forms  of 
displacements  as  given  by  Elquation  (6.14).  Clearly,  in  the  originai  forms  of  equations  (6.2) 
and  (6.3),  the  displacements  show  coupling.  However,  without  this  simplification,  it  would 
be  difficult  to  obtain  the  constant  coefficients  of  the  assumed  displacement  field  in  terms  of 
the  nodal  values  of  displacement,  since  the  difierential  equations  of  Kingsbury  [Kingsbury 
84j  are  set  up  in  terms  of  four  displacements,  while  the  element  has  six  degrees  of  freedom 
per  node.  In  other  words,  in  order  to  obtain  a  finite  element  that  accounts  for  coupling,  we 
would  have  to  impose  coupled  nodal  value  conditions,  which  is  clearly  a  non-trivial  matter. 
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CASE  4: 


«  =  l.OxlO-^'m 
A  =  l.OxlO"'*  /m 

Applied  force  in  the  chord  system,  =  55426.4229645  N 
Computed  Reactions  in  the  Chord  System,  in  N  and  N-m: 


3  Elements 

10  Elements 

F 

Xc 

-  4.4092066  X  10“’ 

-4.4281777x10“’ 

F 

Yc 

-55426.4280001 

-55426.4295408 

F 

Zc 

5.0363461  X  10“  ® 

6.5771379x10“® 

-0.6805012 

-0.6694404 

-6.4613788x10“’ 

-6.7159999x10“’ 

- 121.4721685 

-12.4269992 

Displacements  in  the  Chord  System, 

in  m  and  deg: 

Exact  Solution 

3  Elements 

AX 

C 

5.8975078x10“’ 

5.4112830x10“’ 

<1 

3.9458062x10“’ 

3.9069922x10“’ 

f> 

1.0319861x10“® 

1.4945806x10“® 

^^Xc 

4.2030041  X  10“ 

3.4914457X10“* 

^^Yc 

-4.1399281X10“’ 

2.0471505X10“® 

- 1.2873947 

-1.2035963 

CASE  5: 

K  =  1.0  X  10“'*  /m 
A  =  1.0x10“’  /m 

Applied  force  in  the  chord  system,  =  61441.9525214  N 
Computed  Reactions  in  the  Chord  System,  in  N  and  N>m: 

3  Elements  10  Elements 

F^^  7.4058016  x  10  “  ’  7.4058016  x  lO  “  ’ 


10  Elements 

5.8475781x10“’ 

3.9410363x10“’ 

1.3573915x10“® 

3.3897723x10“* 

2.0838240X10“® 

-1.2799743 
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-61441.9525233 

-61441.9525233 

^Zc 

1.9445444  X  10 

1.9445444x10"* 

-0.7405821 

-0.7405804 

^Yc 

-  7.4658299  X  10"* 

-7.4650474x10"® 

^Zc 

-1.9591481x10"* 

-1.9458587x10"* 

Displacements  in  the  Chord  System, 

in  m  and  deg: 

Exact  Solution 

3  Elements 

10  Elements 

AX 

c 

6.5383257  X  10  "■* 

6.5382582x10-* 

6.5383147x10'* 

ak 

c 

3.6873849x10"* 

3.6873676x10"* 

3.6873678x10'* 

AZ 

c 

1.1442295x  10"® 

1.4901627  X  10"® 

1.4873416x10'® 

5.5909964  x  lO  " 

3.7456808x10'* 

3.7456666x10'* 

<1 

-4.5889408  x  10"® 

2.3109602x10'* 

2.3127627x10'* 

-1.4271180x  10"* 

-1.4286921x10'* 

-1.4287038x10'* 

Comparing  cases  3  and  5,  we  see  that  an  increase  in  the  torsion  has  very  little 
effect  on  the  accuracy  of  finite  element  solution.  On  the  other  hand,  comparing  cases  3  and 
4,  we  see  that  increasing  the  curvature  of  the  gun  tube  reduces  the  accuracy  of  the  solution. 
This  is  consistent  with  the  intuitive  notion  that  when  the  curvature  is  small,  the  effect  of 
torsion  is  to  twist  the  beam  about  its  centroidal  axis  without  changing  its  geometry 
significantly. 


Note  that  even  though  the  displacements  correspond 

to  w,  Ox  and  are  inaccurate,  their  values  are  at  least  an  order  of  magnitude  smaller  than 
the  other  ’’in-plane”  displacements  and  so  the  error  becomes  less  severe.  In  the  next  and 
final  chapter,  we  present  the  conclusions  of  this  investigation. 


CHAPTER  8 


CONCLUSIONS 


In  chapters  2  through  4,  we  presented  an  exact  solution  to  the  problem  of  finding 
the  displacement  of  the  end  of  the  gun  tube  in  terms  of  the  imposed  temperature  change 
that  causes  the  thermal  shroud  to  expand.  In  Chapter  5,  we  performed  a  parametric  study 
of  the  exact  solution  by  plotting  the  components  of  force  and  displacement  against  different 
inititd  curvatures  and  torsions  of  the  gun  tube  for  different  axial  and  flexural  rigidities  of 
the  thermal  shroud. 

A  brief  summary  of  the  conclusions  of  the  parametric  study  follows.  It  was  found 
that  for  a  constant  curvature,  increasing  the  torsion  of  the  gun  tube  had  very  little  effect  on 
the  in-plane  force  and  displacement  components,  while  it  caused  an  increase  in  the  out-of¬ 
plane  force  and  displacement  components.  It  was  observed  that  for  high  and  low  EA, 
displacements  were  minimal.  For  high  El,  the  increased  lateral  resistance  to  bending 
prevented  large  lateral  displacements  of  the  gun  tube.  However,  this  did  not  alleviate  the 
problem  of  relatively  large  angular  displacements.  For  low  EJ,  the  compatibility  condition 
2  reduced  to  condition  1,  which  yielded  fairly  large  displacements.  Thus,  it  appears  that  it 
would  be  most  advantageous  to  have  a  thermal  shroud  with  a  low  EA  for  the  following 
reasons;  Firstly,  the  displacements  are  smallest  in  this  case.  Secondly,  the  forces  are  also 
small.  This  second  reason  makes  the  low  EA  case  preferable  to  the  high  EA  case  since  it 
not  only  reduces  displacements  but  also  the  forces  on  the  connection  between  the  gun  tube 
and  the  shroud.  It  would  be  easier  to  construct  an  axially  flexible  shroud  than  one  which  is 
rigid;  in  addition,  the  gun  tube-shroud  joint  could  be  dengned  for  lesser  loads.  This  would 
lead  to  a  saving  of  material. 


An  alternative  method  of  solution  employing  a  general  space-curved  finite  element 
was  developed  in  Chapter  6.  In  the  preceding  section,  we  compared  the  two  proposed 
solutions.  The  finite  element  solution  yielded  fairly  accurate  results  as  compared  to  the 
exact  solution.  However,  despite  its  relative  inaccuracy,  the  flexibility  of  the  finite  element 
approach  makes  it  highly  attractive.  The  particular  advantages  referred  to  are  the  fact  that 
variable  curvature  and  torsion  can  be  accomodated  within  an  element,  and  also  the 
numerous  possibilities  afforded  in  modeling  a  structure  with  an  assemblage  of  different 
elements.  One  such  example  was  provided  in  Chapter  6,  in  which  the  cross-sections  of 
different  pieces  of  the  gun  tube  were  of  different  sizes.  Thus  it  would  be  possible  to  model  a 
variety  of  highly  complex  structures,  provided  that  the  curvature  and  torsion  of  the 
elements  would  be  small. 


APPENDIX  A 


THE  COMPLEMENTARY  STRAIN  ENERGY  OF  THE  HELICAL  ROD 


The  complete  expression  for  the  streiin  energy  of  deformation,  U,  of  the  helical  rod, 
expressed  in  terms  of  the  XY2  components  of  applied  force  and  moment,  is  as  follows: 
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APPENDIX  B 


THE  STRAIN  ENERGY  OF  DEFORMATION  OF  THE  FINITE  ELEMENT 


The  strain  energy  of  deformation  of  the  space-curved  finite  element  is  a  function  of 
the  nodal  values  of  displacement,  namely, 

U  —  U(u,v.w,u  .v'.0,ii.v.w,u'.v'.<i>  1. 

'  1  ’  l  1  ’  1  ’  1  ’  1’  2’  2  ’  2  ’  2  ’  2  '  ^2  ^ 

The  complete  expression  for  the  strain  energy  is, 

U  =  a3*(l*w2**2+l*wl*w2+l*wl**2)/3,0-a29’'‘((3*l**2*v2p-21»l 

#v2-2*l**2»vlp-9>t‘l*vl)*w2+(2’i'l**2*v2p-0*l*v2-3*l**2*vlp-21*l*vl 
)*wl)/60.0+a30*((l*v2p+6*v2-l*vlp-6*vl)*w2-*‘(-l*v2p+6*v2+l*vlp-6 
*vl)*wl)/12.0+a35*((l*v2p-v2+vl)*w2+(v2-l*vlp-vl)*wl)/l-a26*((3 
*l**2*u2p-21*l*u2-2*l**2*ulp-d*l*ul)*w2+(2*l**2*u2p-9*l*u2-3*l* 
*2*ulp-  21  *l*ul )  *v»l )  /60 . 0+a26*  ( ( I*u2p+6*u2-  l*ulp-6*ul )  *w2+  (  -  l’'‘u2 
p+6*u2+l*ulp-6*ul )  *wl )  /12 .0+a33*  ( (I*u2p-u2+ul )  *v»2+  (u2-l'*ulp-ul ) 

) /l+a49* ( (2*l*p2+l*pl ) *w2+(l*p2+2*l*pl ) *wl) /6 . O+aSS* (w2- wl ) * 
(l*w2+l*wl)/l/2.0+a51*(p2-pl)*(l*w2+l*wl)/l/2.0+a7*(w2-wl)**2/l 
-a31 * (I**2*v2p-6*l*v2-l**2*vlp-6*l*vl ) * (w2-wl ) /1/12 . 0+a36* (v2p- 
Vlp)*(w2*wl)/l+a32*(v2-vl)*(w2-wl)/l-a27*(l**2*u2p-6*l*u2-l**2* 
ulp-6*l*ul)*(w2-wl)/l/12.0+a34*(u2p-ulp)*(w2-wl)/l+a28*(u2-ul)* 
(w2-wl)/l+a50*(l*p2+l*pl)*{w2-wl)/l/2.0+a2*(2*l**3*v2p**2+(-22* 
I**2*v2-3*l**3*vlp-13*l**2*vl)*v2p+78*l*v2**2+(13*l**2*vlp+64*l 
*vl) *v2+2*l**3*vlp**2+22*l**2*vl*vlp+78*l*vl**2)/210.0*al0* (4*1 
J!c*2*v2p**2+{"12*l*v2+4*l**2*vlp+12*l*vl)*v2p+12*v2**2'*(*12*l*vl 
p-24*vl)*v2+4*l**2*vlp**2+12*l*vl*vlp'*-12*vl**2)/l**3*a6*(2*l**2 
*v2p**2+(-3*l*v2-l**2*vlp+3*l*vl)*"2p+18*v2**2*(-3*l*vlp-36*vl) 
*v2+2*l**2*vlp**2+3*l*vl*vlp-*-18*vl**2)/l/16.0-al8*(2*l**2*v2p** 
2+(-18*l*v2-l**2*vlp+3*l*vl)*v2p+18*v2**2+(-3*l*vlp-36*vl)*v2+2 
*l*f2*vlp**2+18*l*vl*vlp*18*vl**2)/l/15.(Hal9*(v2p**2/2.0-vlp** 

2/2 . 0) ♦al 1 ♦ ( (4*l**3*u2p-22*l**2*u2-3*l**3*ulp- 13*l**2*ttl) *v2p+ ( 
-22*l**2*u2p-*'15C*l*tt2+13*l**2*ulp+54*l*ul)*v2+(-3*l**3*ii2p+l3*l 
**2*u2+4*l**3*ulp+22*l**2*ul)*vlp*(-13*l**2*u2p+B4*l*u2+22*l**2 
*ulp+156*l*ul)*vl)/420.0+a20*{(4*l**2*u2p-3*l*u2-l**2*ulp-t3*l*u 
1 )  *v2p+  ( -  3*  I*u2p+36*u2  -3*l*ulp-  36*ul )  *v2+  (  -  l**2*tt2p-  3*  l*u2+4*  1* 
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*2*ulp+3*l*ul)*vlp+(3*l*\i2p-36*u2+3*l*ulp-*-3«*ul)*vl)/l/30.0-a22 
*( (4*l**2*u2p-3*l*u2“l**2*ulp+3*l*ul)*v2p+(”33*l*ii2p+36*u2~3*l* 
ulp-36*ul)*v2+(-l**2*u2p-3*l‘'‘u2'*’4*l**2*ulp+3*l*ul)*vlp+(3*l*u2p 
-36*u2+33*l*ulp+36*ul)*vl)/l/30.0*a37*((4*l**2*u2p-6*l*u2+2*l** 
2»ulp+6*l*ul)»v2p+(-6*l*u2p+12*u2-6*l*ulp-12’''ul)*v2+(2*l**2*u2p 
-6*l*u2+4«l**2*ulp-*'6*l*ul)*vlp+(6*l*u2p-12*u2+6*l*ttlp+12*ul)*vl 
)/l**3-a21*((4*l**2*u2p-33*l*u2-l**2*Tilp+3*l*ul)*v2p+(-3*li'u2p+ 
36*u2-3*l»ulp-36*ul)i'v2*(-l**2*u2p-3*l*u2+4*l**2*ulp+33*l*ul)*v 
lp+(3*l*u2p-36*u2+3*l*ulp+36*ul)*vl)/l/30.0+a23*((l*u2p+2*u2-l* 
ulp-2*ul)*v2p+(2*ulp-2*u2p)*v2+(l*u2p-2*Ti2-l*ulp+2*ul)*vlp+(2*u 
2p-2*ulp)  *vl )  /1/2 . 0+a24*  (  {I*ii2p-2*u2+l*ulp+2*ul )  *v2p+  (2*\i2p-2*\i 
lp)*v2+(-l*u2p+2*u2-l*ulp-2*ttl)*vlp+(2*ulp-2*u2p)*vl)/l/2.0+al3 
*((6*l*u2-l**2*ulp-6*l*ul)*v2p+(-6*l*u2p+30*u2+6*l*ulp+30*ul)*v 
2+(l**2*u2p-6*l*u2+6*l*ul)*vlp+(6*l*u2p-30*u2-6*l*ulp-30*ul)*vl 
) /60 . 0-al2* ( (6*l*u2- l**2*ttlp-6*l*ul ) ♦v2p+ ( -6*l*u2p-30*u2+6*l*ul 
p+30*ul)*v2+(l**2*u2p-6*l*u2+6*l*ul)*vlp-*-(6*l*u2p-30*u2-6*l*ulp 
♦30*ul )  *vl )  /60 . 0-a44*  ( (3*l**2*p2+2*l'<‘*2't‘pl )  *v2p+  ( -21  *l*p2-©*l*p 
l)*v2+(-2*l**2*p2-3*l**2*pl)*vlp+(-9*l*p2-21*l*pl)*vl)/60.0+a45 
*  ( ( l*p2-  l*pl )  *v2p+  (6*p2+6*pl )  *v2+  {l»pl  -  l*p2)  ’''vlp+  ( -6*p2-6*pl )  *v 
l)/12.0+a48*(l*p2’^v2p+(pl-p2)*v2-l*pl*vlp+(p2-pl)*vl)/l-a46*(p2 
-pl)*(l*’''2*v2p-6*l*v2-l**2’'‘vlp-6*l*vl)/l/12.0+al7*(v2**2-vl**2) 
/2.0+a47i‘(p2-pl)*(v2-vl)/l+al*(2*l**3*u2p**2+(-22*l**2*u2-3*l** 
3*ulp-13*l**2*ul)*u2p+78*X*u2**2+(13*l»’*2*ulp+64*l*ul)*u2+2*l** 
3*ulp**2+22*l**2*ul*ulp+78*l*ul**2)/210.0*a9*(4*l**2*u2p**2+(-l 
2*l*u2+4*l**2*ulp+12*l*ul)*u2p+12*ii2**2+(-12*l*ulp-24*ul)''‘tt2+4* 
l**2*ulp*>i<2+12*l*ul*ulp+12''Til’'"*‘2)/l**3+a5*(2*l**2'*u2p’*"*2+(-3*l* 
u2-l**2*ulp+3»l*ul)*u2p'*'18*u2»*2+(-3’'‘l*ulp-36*ul)*u2'*'2’''l**2’''Ulp 
**2+3*l*ul*ulp+18*ul**2)/l/15.0-al5*(2*l**2*u2p**2+(-18*l*u2-l* 
*2*ulp+3*l*ul)*u2p+18*u2*’'‘2+(-3*l*ulp-36*ul)*u2+2*l**2*ulp**2+l 
8*l*Til*ulp+18*ul**2)/l/15.0+al6*(u2p*'»2/2.0-ulp**2/2.0)-a39*((3 
*l**2*p2+2*l**2*pl)*u2p+(-21*l*p2-9*l*pl)*u2+(-2*l**2*p2-3*l**2 
*pl)*ulp+{-9*l*p2-21*l*pl)*ul)/60.(Ha40*((l*p2-l*pl)*u2p+(6*p2+ 
6*pl)*\i2+(l*pl-l*p2)*ulp+(-6*p2-6*pl)*ul)/12.0+a43*(l*p2*u2p'*’(p 
l-p2)*u2-l*pl*ulp+(p2-pl)*ul)/l-a41*(p2-pl)*(l**2*u2p-6*l*u2-l* 
*2*ulp-6*l*ul)/l/12.0+al4*(u2**2-ttl**2)/2.0+a42*(p2-pl)*(u2-ul) 
/l+a4’'‘(l*p2**2+l*pl*p2+l*pl**2)/3.0'*-a8*(p2-pl)**2/l 


where  ul  =  u  ,  ulp  =  u  p  =  ^,  etc.,  and  the  coefBcients  a  ,  a,,  . . . ,  a.,  depend 
on  geometric  and  material  properties  of  the  element,  such  as  the  length  of  the  element,  1, 
the  moments  of  inertia  of  its  cross-section,  ixx,  iyy  and  ixy,  the  curvature,  k  (kappa),  the 
rate  of  change  of  curvature  of  the  element,  kp  (kappa  prinw),  the  torsion,  la  (lambda),  the 
rate  of  change  of  torsion,  lap  (lambda  prime),  and  modulus  of  elasticity  e,  end  finally,  shear 
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modulus  g.  These  coefficients  are  listed  below: 

al  =  (e*(ixx»lap**2+a*k**2)+g*ixx*k**2>t‘la**2)/2.0 

a2  =  iyy*(e*lap**2+g*k**2*la**2)/2.0 

a3  =  iyy*(e*kp**2+g*k**4)/2.0 

a4  =  (ixx*(g*la**2+e*k**2)+g*iyy*la**2)/2.0 

aS  =  (e*ixx*la»*2+g*iyy*k**2)/2.0 

a6  =  (e*iyy*la**2+g*i3cx*k**2)/2.0 

a7  =  e*(iyy*k**2+a)/2.0 

a8  =  g*(iyy+ixx)/2.0 

a9  =  e*iyy/2.0 

alO  =  e*ixx/2.0 

all  =  -ixy*(e*lap**2-*-g*k’»»2*la**2) 
al2  =  -e*ixy*la*lap-g*iyy*k**2*la 
al3  *  g*ixx*k**2*la-a*ixy*la*lap 
al4  =  «*ixx*la*lap+g*ixy*k*»2*la 
al5  =  «*ixy*lap 
al6  =  «*ixy*la 

al7  =  «»iyy*la*lap-g*ixy*k**2*la 

al8  =  -•*lxy*lap 

al9  *  -#*ixy*la 

a20  »  lxy’»(g*k**2-a*la**2) 

a21  s  «*l»c*lap 

a22  =  -a*iyy*lap 
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a23  -  e’^lxx’^la 

a24  =  -e*iyy*la 

a25  =  ixy*(e*kp*lap+g*k*#3*la) 

a26  *  «*ixy*kp*la+g*iyy’t‘k**3 

a27  =  a*ixy*k*lap-a*a*k 

a28  *  e*ixy*k*la 

a29  =  -iyy*(a*kp*lap+g*k**3*la) 

a30  =  g»ixy*k**3-a*iyy*kp’'>la 

a31  =  -e*iyy*k*lap 

a32  =  -e*iyy*k‘i‘la 

a33  =  e*iyy*kp 

a34  *  9*iyy*k 

a35  *  a*ixy*kp 

a36  *  €*ixy*k 

a37  =  «*ixy 

a38  =  a*iyy*k*kp 

a39  =  g*ixy*k*la**2-«*ixx’t‘k*lap 

a40  =  g*iy7*k*la-a*ixx*k*la 

a41  =  g*ixx*k*la 

a42  *  g*ixy*k 

a43  *  -•*ixy*k 

a44  *  a*ixy*k*lap-g*iyy*k*la**2 


a46  »  ixy*(g*k*la+#*k*la) 
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a46  =  -g*ixy*k*la 
a47  =  g*ixx>i‘k 
a48  =  -e*ixx*k 

a49  =  g*iyy*k**2*la-e*ixy*k*kp 
a50  =  -e*ixy*k**2 


a51  =  g*ixy*k**2 


APPENDIX  C 


COMPUTER  IMPLEMENTATION  OF  THE  EXACT  SOLUTION 


C _ 

C  THIS  PROGRAM  COMPUTES  DISPUCEMENTS  AT  THE  END  OF  THE  GUN  TUBE 
C  DUE  TO  THE  FORCE  EXERTED  BY  THE  EXPANSION  OF  THE  THERMAL  SHROUD 

C  FOR  VARYING  INITIAL  CURVATURE  AND  TORSION  OF  THE  GUN  TUBE. 

C _ 

c 

INTEGER  KOUNTl,  K0UNT2.  K0UNT3.  K0UNT4.  lA.  IDGT,  lER 
REAL*8  BLAST,  G.  AC.  AS.  IXX,  JP.  ALPHA.  TMPR,  L.  PI,  CUR,  TOR 
REAL*8  RATIO,  A.  B.  P.  R,  Cl.  SI.  C,  S.  El.  GJ.  EA,  GA 
REAL*8  DENOMl,  DEN0M2,  DEN0M3,  LS.  CVAL,  TVAL 
REAL*8  DISTl,  DIST2.  DISTANCE 

REAL*8  RADI,  RAD2.  RADS.  RA04.  ISHROUD,  EASHROUD.  EISHROUD 
C 

C  DOUBLE  PRECISION  CONSTANTS  ABBREVIATED  FOR  USE  IN  THE  COMPUTATION 
C  OF  FLEXIBILITY  MATRIX  ELEMENTS 

C 

REAL*8  TO.  TR.  FO.  FI.  SI.  SE.  ET.  NI,  TE.  EV,  TW.  FF,  SX.  TF 
REAL*8  TT.  FE.  ST.  NS.  ONE.  ZERO.  START.  FACTOR 
C 

C  ARRAYS  USED: 

C  F  FLEXIBILITY  MATRIX  IN  GLOBAL  COORDINATE  SYSTEM 

C  FBAR  FLEXIBILITY  MATRIX  IN  CHORD  SYSTEM 

C  COMPAT  MATRIX  USED  TO  SOLVE  COMPATIBILITY  EQUATIONS 

C  INVCOMPAT  INVERSE  OF  COMPAT 

C  FORCES  FORCE  VECTOR  IN  COMPATIBILITY  EQUATIONS 

C  RHS  RIGHT  HAND  SIDE  OF  COMPATIBILITY  EQUATIONS 

C  DELTABAR,  DELTA.  DELTATAN 

C  DISPUCEMENTS  IN  CHORD.  GLOBAL  AND  TANGENT  SYSTEMS 

C  FORCEBAR.  FORCE.  FORCETAN 

C  FORCES  IN  CHORD,  GLOBAL  AND  TANGENT  SYSTEMS 

C  CG  CHORD  TO  GLOBAL  TRANSFORMATION  MATRIX 

C  GC  GLOBAL  TO  CHORD  TRANSFORMATION  MATRIX 

C  TG  TANGENT  TO  GLOBAL  TRANSFORMATION  MATRIX 


{ 
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C  GT  GLOBAL  TO  TANGDIT  TRANSFORMATION  MATRIX 

C  TO  TANGENT  TO  CHORD  TRANSFORMATION  MATRIX 

C  TEMP6  TEMPORARY  STORAGE 

C  WKAREA  USED  BY  IMSL  INVERSION  ROUTINE 

C 

REALMS  F(6.6).  FBAR(6.6).  COMPAT(3.3).  INVC0MPAT(3.3) 

REALMS  0ELTABAR(6,1}.  DELTA (6.1).  DELTATAN(6. 1) 

REALMS  FORCEBAR(6.1).  FORCE(6.1).  F0RCETAN(6.1) .  F0RCE3(3.1) 
REAL*8  CG(6.6),  TG(6.6).  CT(6.6).  GC(6.6),  GT(6.6) 

REALMS  TEMP6(6.6).  RHS(3.1).  WKAREACBOO) 

C 

DATA  EUST/4.35007D0/ 

DATA  G/1. 6734509/ 

DATA  RADI.  RAD2.  RAD3.  RAD4/6.0D-2.  l.OD-1,  l.lD-1.  1.45D-1/ 
DATA  ALPHA/1. 206D-5/ 

DATA  TMPR/1.0D2/ 

DATA  L/5.25D0/ 

C 

C  DATA  FOR  CONSTANTS 

C 

DATA  TO.TR.FO.FI.SI.SE/2.ODO.3.ODO.4.ODO.5.ODO.6.ODO.7.ODO/ 

DATA  ET.ni.TE.EV.TW.FF/8.0D0.9.0D0.10.0D0.11.0D0.12.0D0.16.0D0/ 
DATA  SX.TF.TT.FE.ST/16.0D0.24,OD0.33.0D0.48.0D0.72.0DO/ 

DATA  NS . ONE . ZERO . START/96 . ODO . 1 . ODO . 0 . ODO . 1 . 0D*4/ 

C 

C  OPEN  OUTPUT  FILES 

C 

0PEH(UNIT=12 .FILE* ’dtlxbar.dat ‘ .STATUS* ’NEW* ) 

OPEN (UNIT*13 .FILE* ’dalybar .dat ’ .STATUS* ’NEW’ ) 

OPEN (UNIT* 14 .FILE* ’dalzbar.dat ’ .STATUS* ’NEW’ ) 

OPEN (UNIT* 1 5 .FILE* ’ thaxbar . dat ’ . STATUS* ’ NEW ’ ) 

OPEN (UNIT*16 .FILE* ’ thaybar . dat ’ . STATUS* ’ NEW ’ ) 

OPEN (UNIT*17 .FILE* ’ thazbar . dat ’ . STATUS* ’NEW ’ ) 

OPEN (UNIT*34 .FILE* ’forcaxbar.dat ’ . STATUS*’ NEW ’ ) 

0PEN(UNIT*18 .FILE* ’lorcaybar .dat ’ .STATUS* ’NEW ’ ) 

OPEN (UNIT*37 . F ILE* ’ 1 orcazbar . dat ’ .STATUS* ’ NEW ’ ) 

OPEN (UNIT*33 .FILE* ’dlatanca.dat’ .STATUS* ’NEW’ ) 

C 

C  ASSIGN  VALUES  OF  PI  AND  FACTOR  (FOR  MULTIPLYING  CURVATURE) 

C 

PI  *  OATAN(ONE)*FO 

FACTOR  «  OSqRT(OSqRT(OSqRT(DSQRT(TE)))) 

C 

C  COMPUTE  SECTION  PROPERTIES  OF  GUN  TUBE  AND  THERMAL  SHROUD 
C 
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C  GUN  TUBE: 

C 

AC  =  PI*(RAD2**2-RAD1**2) 

IXX  »  PI*(RAD2**4-RAD1**4)/F0 
JP  =  TO* IXX 
El  »  ELAST*IXX 
GJ  =  G*JP 
EA  =  ELAST*AC 
GA  =  G*AC 
C 

C  THERMAL  SHROUD; 

C 

AS  =  PI*(RAD4**2-RAD3**2) 

ISHROUD  =  PI*(RAD4**4-RAD3**4)/F0 
EISHROUD  =  ELAST* ISHROUD 
EASHROUD  =  ELAST*AS 
C 

C  INITIALIZE  TOR  AND  TVAL 

C 

TOR  =  ZERO 
TVAL  =  ZERO 
C 

C  LOOP  TO  COMPUTE  DISPLACEMENTS  FOR  VARYING  CURVATURES  AND  TORSIONS 
C 

DO  100  K0UNT3  *  1 .  50 
CUR  s  START 
CVAL  =  ONE 

IF  (K0UNT3.Eq.2)  THEN 
TOR  =  START 
TVAL  =  ONE 
END  IF 

DO  200  K0UNT4  *  1 .  49 
C 

C  COMPUTE  HELIX  GEOMETRY  FROM  CUR  AND  TOR  AND  INITIALIZE  ARRAYS 
C 

RATIO  »  TOR/CUR 
A  =  L/DSQRT(0NE+RATI0**2) 

B  =  DSqRT(L**2-A**2) 

Cl  »  A/L 
SI  »  B/L 
R  «  C1**2/CUR 
P  •  A/R 
C  «  DCOS(P) 

S  «  DSIN(P) 

LS  »  DSQRT((R*C-R)**2*(R*S)**2*B**2) 


130 


C 

C  INITIALIZE  TRANSFORMATION  MATRICES  TO  ZEROES 
C 

DO  101  KOUNTl  »  1.  6 
DO  101  K0UNT2  =1.6 

CG (KOUNTl.  K0UNT2)  =  ZERO 
TG (KOUNTl.  K0UNT2)  *  ZERO 
CT(K0UNT1.  K0UNT2)  =  ZERO 
GC (KOUNTl.  K0UNT2)  =  ZERO 
GT(K0UNT1,  K0UNT2)  *  ZERO 

101  CONTINUE 
C 

C  INITIALIZE  FORCE  VECTORS  IN  ALL  3  COORDINATE  SYSTEMS  TO  ZEROES 
C 

DO  102  KOUNTl  =1.6 

FORCEBAR (KOUNTl .  1)  =  ZERO 
F0RCE(K0UNT1 .  1)  »  ZERO 
FORCETAN (KOUNTl.  1)  =  ZERO 

102  CONTINUE 
C 

C  COMPUTE  ALL  6  TRANSFORMATION  MATRICES 

C  FIRST.  ASSIGN  ELEMENTS  OF  CG  (CHORD  TO  GLOBAL  TRANSFORMATION  MATRIX) 
C 

DENOMl  =  DSqRT(TO*(ONE-C)) 

DEN0M2  =  0SqRT(T0*(0HE-C)+(P*B/A)**2) 

0EN0M3  =  DENOMl *DEN0M2 
CG(l.l)  =  S/DENOMl 
CG(1.2)  =  (C-0NE)/DEN0M2 
CG(1.3)  =  (0NE-C)*(P*B/A)/DEN0M3 
CG(2.1)  =  (ONE-C) /DENOMl 
CG(2.2)  =  S/DEN0M2 
CG(2.3)  =  -(S*P*B/A)/DEJI0M3 
CG(3.1)  =  ZERO 
CG(3.2)  *  (P*B/A)/DEN0M2 
CG(3.3)  =  TO* (ONE-C) /DEN0M3 
DO  103  KOUNTl  =1.3 
DO  103  K0UNT2  =1.3 

CG(K0UNTl+3.  K0UNT2+3)  =  CG(K0DNT1.  K0DNT2) 

103  CONTINUE 
C 

C  ASSIGN  ELEMENTS  OF  TG  (TANGENT  TO  GLOBAL  TRANSFORMATION  MATRIX) 

C 

TG(l.l)  =  -C 
TG(1.2)  =  S1*S 
TG(1.3)  =  -C1*S 
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TG(2,1)  =  -S 
TG(2,2)  =  -S1*C 
TG(2,3)  =  C1*C 
TG(3.1)  =  ZERO 
TG(3.2)  =  Cl 
TG(3,3)  *  SI 
DO  104  KOUNTl  =1.3 
DO  104  K0UNT2  =1.3 

TG(K0UNTl+3.  K0UNT2+3)  =  TG(K0UMT1.  K0UNT2) 

104  CONTIIIUE 
C 

C  INVERT  CG  TO  GET  GC  (GLOBAL  TO  CHORD  TRANSFORMATION  MATRIX) 

C  AND  TG  TO  GET  GT  (GLOBAL  TO  TANGENT  TRANSFORMATION  MATRIX) 

C 

DO  106  KOUNTl  =1.6 
DO  106  K0UNT2  =1.6 

TEMP6 (KOUNTl .  K0UNT2)  =  CG(KOUNTl.  K0UNT2) 

106  CONTINUE 
IDGT  =  0 
lA  =  6 
KOUNTl  =  6 

CALL  LINVIF  (TEMP6 .  KOUNTl .  I A .  GC .  IDGT.,  WKAREA ,  lER) 

DO  106  KOUNTl  =1.6 
DO  106  K0UNT2  =1.6 

TEMP6(K0UNT1 .  K0UNT2)  »  TG(K0UNT1.  K0UNT2) 

106  CONTINUE 
IDGT  =  0 
lA  =  6 
KOUNTl  =  6 

CALL  LINVIF (TEMP6 .KOUNTl , lA . GT . IDGT . WKAREA . lER) 

C 

C  COMPUTE  CT  (CHORD  TO  TANGENT  TRANSFORMATION  MATRIX)  »  GT  *  CG 
C 

KOUNTl  *  6 

CALL  MULMATX(GT.  CG.  CT.  KOUNTl.  KOUNTl.  KOUNTl) 

C 

C  ASSIGN  ELEMENTS  OF  F  MATRIX 
C 

F(l,l)  ■  L**3/(TW*EI)*(SI*C1**4*(TR*S*C“F0*S+T0*P*S**2+P) 

1  ♦S1**4*(TR*P*T0*P**3'TR*S*C)+SI*C1**2*S1**2*(F0*S- 

2  TR>»S*C-P)  ••■Sl**2*  (-TR*P*I0*P**3+TR*S*C)  )  /P**3 

3  -L**3/(TF*GJ)*C1**2*S1'**2/P**3*(-FF*P-TW*P*S**2 

4  “T0*P**3+FE*S-TT*S*C) 

6  ♦L/(T0*EA)*C1**2*(P-S*C)/P  +  L/(F0*GA)*(Cl**2*(P+ 
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6  S*C)/P+T0*S1**2) 

F(l,2)  *  L**3/(F0*EI)*(T0*C1**4*(TR*S**2-T0*P’»S*C4T0>*‘C-T0)  + 

1  S1**4*(P**2-S**2)+T0*C1**2*S1**2*(-TR*S**2+T0*P*S 

2  -T0*C+T0)+S1**2*(S**2"P**2))/P**3 

3  -L**3/(ET*GJ)*C1**2*S1**2*(ET“ET*C-EV*S**2+F0*P*S*C 

4  +F0*P*S'P**2)/Pt‘*3 

6  -L/(T0*EA)*Cl**2*S*=t‘2/P  L/(F0*GA)*C1**2*S**2/P 

?(i.3)  =  L**3/(F0*El)*(Cl*Sl*P*(S-P*C)+T0*Cl‘i‘*3*Si*(-T0+T0*C 

1  -T0*S**2+TR*P’'‘S)+C1*S1**3*(F0*F0*C-P*S-P**2*C))/P**3 

2  -L**3/(ET*GJ)’t‘Cl**3'''Sl*(P’'‘*2*C+SE'*P*S-FQ*S**2+ET*C* 

3  ET)/P**3-L/EA*Cl>fSl*(0NE-C)/P+L/(T0*GA)*Cl*Sl*(0ME- 

4  C)/P 

F(l,4)  *  L**2/(F0*EI)*(S1*(P-S*C)+T0*C1**2*S1*(P+S*C*T0*S)+ 

1  S1**3*(S*C-P))/P**2 

2  -L**2/{ET*GJ)*C1**2*S1*(S*C-F0*S+TR*P)/P**2 

F(1.5)  «  L**2/(F0*EI)*(Sl*(P**2-S**2)-«-T0*Cl**2*Sl*S**2+Sl**3 

1  *(p*f2+S**2))/P**2 

2  •L**2/(ET*GJ) ♦01**2*81* (S**2-P**2) /P**2 

F(l,6)  a  L**2/EI*(C1**3*(0ME-C-P*S)+C1*S1**2*(C-0NE))/P**2 
1  -L**2/(T0*GJ)*C1*S1**2*{P*S+T0*C-T0)/P**2 

F(2.2)  a  L**3/(TW*EI)*(SI*C1**4*(P-TR*S*C+T0*P*C**2)+SI* 

1  C1**2*S1**2*(P+TR*S*C-F0*P*C)+S1**4*(-TR*P+T0*P**3+ 

2  TR*S*C)  +81  **2*  (TR*P+T0*Paa3*TR*S*C)  )  /P**3 

3  -L**3/ {TF*GJ) *C1**2*81**2* ( *NIaP-TO*P**3-TF*P*C+TW* 

4  p*S**2+TT*8*C)/P**3 

6  +L/(T0*EA)*C1**2*(P+8*C)/P  +  L/(F0*GA) *(01**2* (P-8 

6  *C)/P+T0*81**2) 

F(2,3)  a  L**3/(F0*EI)*(Cl*81*(P*C-8-P**2*S)+T0*Cl**3*81*(8+ 

1  T0*8*C-TR*P*C)+C1*81**3*(-TR*S-P*C-P**2*S+F0*P))/ 

2  P**3 

3  -L**3/(ET*GJ)*C1**3*81*(FI*S-FI*P*C+F0*S*C+P**2*8 

4  -F0*P)/P**3  ♦  L/EA*C1*S1*S/P  -  L/(T0*GA)*C1*S1*S/P 

F(2.4)  a  L**2/(F0*EI)*(81*(-S**2-P**2)+T0*C1**2*S1*(-C**2+I0 

1  ♦C-0NE)+81**3*(S**2-P*a2))/P**2 

2  -L**2/(ET*GJ)*C1**2*S1*(-F0+F0*C+S**2+P**2)/P**2 

F(2,6)  a  L**2/(FO*EI)*(-81+TO*C1**2*S1+S1**3)*(P-8*C)/P**2 
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1  -L**2/(ET*GJ)*C1**2*S1*(P-S*C)/P**2 

F(2.6)  *  L**2/EI*(C1**3*(P*C-S)+C1*S1**2*(S-P))/P**2 
1  -L**2/ (TO*GJ) *C1*S1**2*  C~P*C-P*TO*S) /P**2 

F(3.3)  *  L**3/(T0*EI)*(C1**2*(P-S*C)+C1**2*S1**2*(TR*P+S*C- 

1  F0*S))/P**3 

2  -L**3/(F0*GJ)*C1**4*(F0*S-TR*P-S*C)/P**3 

3  +L/EA*S1**2  +  L/(T0*GA)*C1**2 

F(3.4)  =  L**2/(T0*EI)*(C1*P*S+C1*S1**2*(T0*C-T0+P*S))/P**2 
1  -L**2/ (FO»GJ) *01**3* (-TO^C+TO-P+S) /P**2 

F(3.5)  =  L**2/(T0*EI)*(C1+C1*S1**2)*(S-P*C)/P**2 
1  -L**2/(F0*GJ)*C1**3*(P*C-S)/P**2 

F(3.6)  =  L**2*(0NE/EI-0NE/(T0*GJ))*C1**2*S1*(S-P)/P**2 

F(4.4)  =  L/(T0*EI*P)*(S*C*C1**2+P*(0NE+S1**2)) 

1  -L*C1**2/(F0*GJ*P)*(S*C-P) 

F(4.5)  =  a/(T0*EI)-L/(F0*GJ))*{Cl**2*S**2/P) 

F(4.6)  =  a/EI-L/(TO*GJ))*(Cl*Sl*(ONE-C)/P) 

F(5.5)  *  L/(T0*EI*P)*(-C1**2*S*C+P*(0NE+S1**2)) 

1  -L/(F0*GJ*P)*(C1**2*(-S*C-P)) 

F(5.6)  =  (L/(T0*GJ)-L/EI)*(C1*S1*S/P) 

F(6.6)  =  L*C1**2/EI+L*S1**2/(T0*GJ) 

DO  108  KOUNTl  =2.6 
DO  108  K0UNT2  »  1,  KOUNTl -1 

FCKOUNTl,  K0UNT2)  =  F(K0UNT2.  KOUNTl) 

108  CONTINUE 
C 

C  PRE  *  POSTMULTIPLY  ORIGINAL  F  WITH  GO  ft  CG  RESPECTIVELY  TO 
C  OBTAIN  THE  TRANSFORMED  FLEXIBILITY  MATRIX.  FBAR.  IN  THE  CHORD  SYSTEM 
C 

KOUNTl  «  6 

CALL  MULMATX(GC.F.TEMP6,K0UNT1.K0UNT1 .KOUNTl) 

KOUNTl  «  6 

CALL  MULMATX(TEMP6.CG.FBAR.K0UNT1.K0UNT1 .KOUNTl) 


C 
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C  APPLY  COMPATIBILITY  OF  DISPUCEMENTS  TO  FIND  FORCES  EXERTED  ON 
C  GUIl  TUBE  DUE  TO  EXPANSION  OF  SHROUD 
C 

C  COMPATIBILITY  IS  EXPRESSED  AS  A  SET  OF  3  SIMULTANEOUS  EQUATIONS 

C  IN  THE  3  FORCE  COMPONENTS.  FXBAR.  FYBAR,  FZBAR 

C 

COMPATd.l)  »  FBAR(l,l)+LS**3/(TR*EISHROUD) 

C0MPAT(1.2)  *  FBAR(1.2) 

C0MPAT(1.3)  «  FBAR(1.3) 

C0MPAT(2.1)  *  FBAR(2.1) 

C0MPAT(2.2)  «  FBAR(2.2)+LS/(EASHR0UD) 

C0MPAT(2.3)  »  FBAR(2,3) 

C0MPAT(3.1)  *  FBAR(3.1) 

C0MPAT(3.2)  =  FBAR(3.2) 

C0MPAT(3.3)  =  FBAR(3.3)+LS**3/(TR*EISHR0UD) 

RHS(l.l)  *  ZERO 
RHS(2.1)  =  LS*ALPHA*TMPR 
RHS(3.1)  »  ZERO 

IDGT  =  0 
lA  =  3 
KOUNTl  »  3 

CALL  LINVlFCCOMPAT.KOOTfTl .lA.INVCOMPAT.IDGT.WKAREA.IER) 

KOUNTl  =  3 
K0UNT2  =  1 

CALL  MULMATXdNVCOMPAT .RHS ,F0RCE3. KOUNTl .KOUNTl .K0UNT2) 

DO  110  KOUNTl  »  1.  3 

FORCEBAR (KOUNTl. 1)  »  F0RCE3 (KOUNTl. 1) 

110  CONTINUE 
C 

C  OBTAIN  DISPLACEMENTS  VECTOR  AS  DELTABAR  -  FBAR  *  FORCEBAR 
C 

KOUNTl  *  6 
K0UNT2  »  1 

CALL  MULMATX(FBAR.F0RCEBAR.DELTABAR.K0UNT1 .KOUNTl .K0UNT2) 

C 

C  CONVERT  DISPUCEMENTS  AND  FORCES  INTO  DIFFERENT  SYSTEMS 
C 

KOUNTl  «  6 
K00NT2  »  1 

CALL  MULMATX(C6.F0RCEBAR.F0RCE.K0UNT1 .KOUNTl .K0UNT2) 

CALL  MULMATX(CT. FORCEBAR. FORCETAN.KOUNTl .KOUNTl .K0UNT2) 
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CALL  MULMATX(CG .DELTABAR, DELTA. KOUNTl .KOUNTl .K0UNT2) 
CALL  MULMATX (CT .DELTABAR .DELTATAN .KOUNTl .KOUNTl .K0UNT2) 
C 

C  COirVERT  ANGULAR  DISPLACEMENTS  INTO  DEGREES 
C 

DO  109  KOUNTl  «  4.  6 

DELTABAR (KOUNTl. 1)  =  DELTABAR(K0UNT1 .1)*1 .8D2/PI 
DELTA (KOUNTl. 1)  =  DELTA(K0UNT1.1)*1.8D2/PI 
DELTATAN(K0UNT1 . 1 )  =  DELTATAN(K0UNT1 . 1 ) *1 . 8D2/PI 
109  COKTIlfUE 

DISTANCE  =  DSQRT(DELTABAR(1,1)**2  ♦  DELTABAR(2. 1)**2  + 

1  DELTABAR (3. 1)**2) 

DISTl  =  DSQRT(DELTA(1.1)**2  +  DELTA(2.1)**2  + 

1  DELTA(3.1)**2) 

DIST2  =  DSQRT (DELTATAN (1.1)**2  ♦  DELTATAN (2. 1)**2  + 

1  DELTATAN(3.1)**2) 

C 

C  'JlRITE  THE  RESULTS  IN  THE  OUTPUT  FILES  OPENED 
C 

WRITE (12.*)  CVAL.  TVAL,  DELTABARd .1) 

WRITE (13.*)  CVAL.  TVAL.  DELTABAR(2.1) 

WRITE (14.*)  CVAL.  TVAL.  DELTABAR(3.1) 

WRITEdS.*)  CVAL.  TVAL.  DELTABAR(4.1) 

WRITEde.*)  CVAL.  TVAL.  DELTABAR(6. 1) 

WRITEd?.*)  CVAL.  TVAL.  DELTABAR(6. 1) 

WRITE (34.*)  CVAL.  TVAL.  FORCEBARd.l) 

V®ITEd8.*)  CVAL.  TVAL.  F0RCEBAR(2.1) 

WRITE(37.*)  CVAL,  TVAL.  F0RCEBAR(3.1) 

WRITE(33.*)  CVAL.  TVAL.  DISTANCE 
C 

C  START  OVER  WITH  NEW  CUR  AND  TOR  VALUES 
C 

CUR  *  CUR*FACTOR 
CVAL  =  CVAL  +  6.2SD-2 
200  CONTINUE 

TOR  *  TOR*FACTOR 
TVAL  =  TVAL  *  6.26D-2 
100  CONTINUE 
STOP 
END 
C 

C  MULMATX  SUBROUTINE  BEGINS 
C 

SUBROUTINE  MULMATX(MAT1 ,MAT2.MAT3,L1 ,L2,L3) 

INTEGER  LI.  L2,  L3.  KNTl.  KNT2.  KNT3 
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REAL*8  MAT1(L1.L2).  MAT2(L2.L3),  MAT3(L1.L3) 

DO  111  KNTl  =  1.  LI 
DO  111  KNT2  *  1.  L3 

MAT3(KNT1.KNT2)  »  O.ODO 
DO  112  KNT3  =  1.  L2 

MAT3(KNT1.KNT2)  =  MAT3(KNT1.KNT2)+MAT1(KNT1,KNT3)* 
1  MAT2(KNT3.KNT2) 

112  CONTINUE 
111  CONTUSE 
RETURN 
END 


C 

C 

C 

C, 


END  OF  PROGRAM 


APPENDIX  D 


COMPUTER  IMPLEMENTATION  OF  THE  FINITE  ELEMENT  METHOD 


C _ 

c 

C  DECLARATIONS 

C 

C _ 

c 

c 

INTEGER  I.  lA.  IDGT,  lER.  J.  NELEM.  NI.  NJ.  NK 


C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


REALMS  K.  KP.  LA.  LAP.  A.  IXX.  lYY.  IXY.  E.  G 
REALMS  AA.  BB.  RATIO.  P.  PP.  R.  S.  C.  Si.  Cl.  L.  LTUBE 
REALkB  DEMOMl.  DEN0M2.  DEN0M3.  PI.  RINNER.  ROUTER 

REALMS  XO.  XI.  X2.  X3.  X4.  X5.  X6,  X7,  XO.  XIO.  Xll 
REAL*8  X12.  X13.  X15.  X20.  X30.  X35.  X60.  X70,  X105 
REALMS  X140.  X180.  X210.  X420 

ARRAYS  USED: 


CG 

GC 

TG 

TC 

TCINV 

TRANS 

TRANSINV 

ELSTIF 

GLSTIF 

CHDSTF 

REDFLEX 

FKNOWN 

OELCHD 


CHORD  TO  GLOBAL  TRANSFORMATION  MATRIX 

GLOBAL  TO  CHORD  TRANSFORMATION  MATRIX 

TANGENT  TO  GLOBAL  TRANSFORMATION  MATRIX 

TANGENT  TO  CHORD  TRANSFORMATION  MATRIX 

INVERSE  OF  TC 

GLOBAL  TC  MATRIX 

INVERSE  OF  GLOBAL  TC  MATRIX 

ELEMENT  STIFFNESS  MATRIX 

GLOBAL  STIFFNESS  MATRIX 

GLOBAL  STIFFNESS  MATRIX  IN  CHORD  SYSTEM 

REDUCED  FLEXIBILITY  MATRIX 

GLOBAL  FORCE  VECTOR  IN  CHORD  SYSTEM 

GLOBAL  OISPUCEMENT  VECTOR  IN  CHORD  SYSTEM 


REAL*8  CG(3,3).  GC(3.3).  TG(3,3).  TC(3,3).  TCINV(3.3) 
REAL’^8  TRANS(24.24).  TRANSINV(24.34) .  ELSTIF(12.12) 
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138 


C 

C 

C 


C 


C, 

C 

C 

C 

C 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 


REALMS  GLSTIF(24.24).  CHDSTF(24 .24) .  REDFLEX(18.18) 

REAL*8  WKAREA(IOO).  TEMP61(6.1).  TEMP18(18.18) ,  TEMP24(24.24) 
REAL*8  FKN0W1I(18.1) .  DELCHD(18. 1) .  CHECK(18,18) 

REAL*8  KBETAALPHA(6.18) 

DOUBLE  PRECISION  CONSTANTS  USED  IN  STIFFNESS  MATRIX  DEFINITION 

DATA  X0.X1,X2.X3.X4/0.0D0.1.0D0,2.0D0,3.0D0.4.0D0/ 

DATA  X5.X6.X7.X®,X10/5.0D0.6.0D0.7.0D0.9.0D0.1.0D1/ 

DATA  X11.X12.X13.X16.X20/1.1D1.1.2D1.1.3D1.1.BD1,2.0D1/ 

DATA  X30 , X35 . X60 . X70 . X106/3 . ODl . 3 . 5D1 .6 . ODl . 7 . ODl . 1 . 05D2/ 

DATA  X 140 .  X180 .  X210 .  X420/ 1 .  4D2 . 1 .  8D2 . 2 . 11)2 . 4 . 2D2/ 

OPEI^  (UNIT=  1 1 .  FILE=  ’  f  €mdata  ’ .  STATUS=  *  OLD  ’ ) 

OPEN (UNIT= 1 2 . FILE= ’ f amoutput ’ .STATUS=  *  NEW ’ ) 


ASSIGN  VALUE  OF  PI  USING  THE  FACT  THAT  TAN (PI/4)  =  1 
PI  =  DATA1I(X1)*X4 


(I)  READ  MATERIAL  AND  GEOMETRIC  PROPERTIES  OF  EACH  ELEMENT 

(II)  COMPUTE  ELEMENT  STIFFNESS  MATRIX 

(III)  ADD  THE  CONTRIBUTION  OF  THE  ELEMENT  TO  THE  GLOBAL 
STIFFNESS  MATRIX 

THIS  ASSUMES  CONNECTIVITY  INFORMATION.  I.E..  THE  ELEMENTS  ARE 
CONNECTED  TO  EACH  OTHER  END  TO  END  IN  THE  ORDER  GIVEN 


READ  GEOMETRICAL  AND  MATERIAL  PROPERTIES 
READdl.*)  NELEM,  K.  KP.  LA.  UP.  E.  G.  LTUBE 
COMPUTE  GEOMETRICAL  PARAMETERS  OF  THE  HELIX 
RATIO  -  U/K 

AA  «  LTUBE/DSqRT(Xl  +  RATI0**2) 

BB  >  DSqRT(LTUBE<f*2  •  AA**2) 

Cl  «  AA/LTUBE 
SI  «  BB/LTUBE 
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R  =  Cl**2/K 
P  =  AA/R 
C  =  DCOS(P) 

S  =  DSIN(P) 

C 

C  ASSIGN  VALUES  TO  CG  (CHORD  TO  GLOBAL  TRANSFORMATION  MATRIX) 

C 

DENOMl  =  DSqRT(X2*(Xl-C)) 

DEN0M2  =  DSQRT(X2*(X1-C)  +  (P*BB/AA)**2) 

DEN0M3  =  DENOMl *DEN0M2 

CG(l.l)  =  S/DENOMl 

CG(1.2)  =  (C-X1)/DEN0M2 

CG(1.3)  =  (X1-C)*(P<<BB/AA)/DEN0M3 

CG(2.1)  =  (XI -C) /DENOMl 

CG(2.2)  =  S/DEN0M2 

CG(2.3)  =  -(S*P*BB/AA)/DEN0M3 

CG(3.1)  =  XO 

CG(3.2)  =  (P*BB/AA)/DEN0M2 
CG(3.3)  =  X2»(X1-C)/DEM0M3 
C 

C  IirVERT  CG  TO  OBTAIN  GC  (GLOBAL  TO  CHORD  TRANSFORMATION  MATRIX) 

C 

IDGT  =  0 
lA  =  3 
NI  »  3 

CALL  LINV1F(CG.  NI .  lA.  GC.  IDGT.  WKAREA.  lER) 

C 

C  INITIALIZE  ELEMENT  AITO  GLOBAL  STIFFNESS  MATRICES.  TRANSFORMATION 

C  MATRICES.  GLOBAL  FORCE  VECTOR.  CONNECTIVES  I  AND  J.  AND  ANGLE  PP 

C  SUBTENDED  BY  ARC  AT  GLOBAL  ORIGIN 
C 

DO  100  NI  =  1.  12 
DO  100  NJ  =  1.  12 
ELSTIF(NI.NJ)  =  XO 

100  CONTINUE 

DO  101  NI  =  1.  6*(NELEM+1) 

DO  101  NJ  =  1.  6*(NELEM+1) 

GLSTIF(HI.NJ)  =  XO 
TRANS(NI,NJ)  =  XO 
TRANSINV(NI.NJ)  =  XO 

101  CONTINUE 

DO  102  NI  >  1.  6*NELEM 
FKNOWN(NI.l)  =  XO 

102  CONTINUE 
1  =  0 


i 
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J  *  0 
PP  *  XO 
C 

C  READ  APPLIED  FORCES  AIR)  MOMENTS 
C 

READCll,*)  (FKN0WN(6*(NELEM-l)+NI.l).  NI  -  1,  6) 

C 

C  FOR  EACH  ELEMENT.  DO  THE  FOLLOWING: 

C 

DO  200  NK  =  1.  NELEM 

READdl.*)  RINNER,  ROUTER.  L 

A  =  PI*(R0UTER**2-RINNER**2) 

IXX  =  PI*(R0UTER**4-RINNER**4)/X4 
lYY  =  IXX 
IXY  =  XO 

ELSTIF(l.l)  =  X13*L*(E*(IXX*LAP*#2-*-A*K**2)+CfIXX*K**2*LA**2)/X35 

1  +X12’i‘((E*IXX^LA**2+G*IYY*K**2)/X2-E*IXY*UP)/(X5*L)-E*IXX*LA 

2  *LAP-G-^IXY>!‘K**2*LA+X12*E*IYY/L*=*'3 

ELSTIF(1.2)  =  (-X13)*IXY*L*(E*LAP**2+G*K**2*U**2)/X36-(-X2*E* 

1  IXY*LA*LAP-G*IYy’i‘K**2*LA+G*IXX*K**2*LA)/X2+(-X6)*(*E*IYY*LAP+ 

2  E*IXX*LAP- IXY* (G*K**2-E*U*»2) ) / (XB*L) ♦X12*E*IXY/L**3 

ELSTIF(1.3)  =  X7*IXY*L*(E*KP*LAP+G*K**3*LA)/X20-(E*IXY*K*UP+E* 

1  IXY*KP*LA+G*iyY*K**3-A*E*K)/X2+(E*IXY*K*U-E*IYY*KP)/L 

ELSTIF(1.4)  =  X11*IXY*L**2*(E*LAP**2+G*K**2*LA**2)/X210-(E*IYY* 

1  LAP-X11*E*IXX*LAP+IXY*(G*K**2-E*LA**2))/X10-L*(G*IYY*K**2*LA 

2  +G*IXX*K**2*U)/X10-(E*IYY*U+E*IXX*U)/L-X6*E*IXY/L**2 

ELSTIFd.S)  =  Xll*L**2*(E*(IXX*LAP**2t-A*K**2)+C*IXX*K**2*U**2) 

1  /X210+( (E*IXX*LA**2+G*IYY*K**2)/X2-X6*E*IXY*LAP)/X5+X6*E*IYY 

2  /L**2 

ELSTIFd.6)  *  X7*L*(G*IXY*K*U**2*E*IXX*K*LAP)/X20-(G*1YY*K*U* 

1  G*IXX*K*U-E*IXX*K*U)/X2+(G*IXy*K+E*IXY*K)/L 

ELSTIFd.7)  »  X9*L*(E*(IXX*LAP**2+A*K**2)+G*IXX*K**2*U**2)/X7<H 

1  X12*(E*IXY*LAP-(E*IXX*U**2+G*IYY*K**2)/X2)/(X6*L)-X12*E*ITY/ 

2  L**3 


ELSTIFd.S)  «  (-X9)*IXY*L*CE*LAP**2+G*K**2*U**2)/X70*X6*(-E*IYY 
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1  *LAP+E*  m*LAP-  IXY*  (G*K**2-E*LA**2) )  /  (X5*L) + (G*IYY*K**2*LA+G* 

2  IXX#K**2*LA)/X2-X12*E*IXy/L**3 

ELSTIFCl.G)  =  X3*IXY*L*(E*KP*LAP>>G*K**3*LA)/X20+(E*1XY*K*LAP-E* 

1  IXY*KP*LA-G*IYY*K**3-A*E*K)/X2+(E*iyY*KP-E*IXY*K*U)/L 

ELSTIFCl.lO)  =  (-X13)*IXY*L**2*(E*LAP**2+G*K**2*U**2)/X420-(E* 

1  IYY*LAP-E*IXX*UP+IXY*(G*K**2-E*LA**2))/X10-L*(-GfIYY*K**2*U 

2  -G*IXX*K**2*LA)/X10-(-E*IYY*U-E*IXX*LA)/L-X6*E*IXY/L**2 

ELSTIFCl.ll)  =  (-X13)*L**2*(E*(IXX*LAP**2*A*K**2)*G*IXX*K**2* 

1  LA**2)/X420+((E*IXX*LAf'*2+G*IYY*K**2)/X2-E*IXY*LAP)/X5+X6*E 

2  *1YY/L**2 

ELSTIF(1.12)  =  X3*L*(G*rXY*K*LA**2-E*IXX>»K*LAP)/X20+(-G*IYY*K*U 
1  +G* IXX*K*LA+E* IXXkK^LA) /X2+ ( -G*IXY*K-E* IXY*K) /L 

ELSTIF{2.2)  =  X13*IYY*L*(E*LAP**2+G*K**2*LA**2)/X35+X12*(E*IXY* 

1  LAP+(E*IYY*LA**2+G*IXX*K**2)/X2)/(XB*L)-E*IYY*LA*LAP+G*IXY* 

2  K**2*LA+X12*E*IXX/L**3 

ELSTIF(2.3)  =  (-X7)*IYY*L*(E*KP*UP+G*K**3*LA)/X20-(-E*IYY*K*UP 
1  -E*IYY*KP*U+G*IXY*K»*3)/X2-^(-E*IYY*K*U-E*IXY*KP)/L 

ELSTIF(2.4)  =  (-X11)*IYY*L**2*{E*UP**2+G*K**2*LA**2)/X210-(X6*E* 
1  IXY*LAP+(E*1YY*LA**2+G*IXX*K**2)/X2)/X6-X6*E*IXX/L**2 

ELSTIF (2,5)  =  (-X11)* IXY*L**2* (E»LAP**2+G*K**2*U**2) /X210+ (XI 1 *E 

1  *IYY*UP-E*IXX*LAP+IXY*(G''‘K**2-E*LA**2))/X10+L*(-G*IYY*K**2*U 

2  -G*IXX*K**2*U)/X10+(-E*IYY*U-E*IXX*LA)/L+X6*E*IXY/L**2 

ELSTIF(2.6)  =  X7*L*(E*IXY*K*LAP-G*IYY*K*U**2)/X20-(IXY*(G*K*U+ 

1  E*K*U)-G*IXY*K*LA)/X2+(G*IXX*K+E*IXX*K)/L 

ELSTIF(2.7)  «  (-X0)*IXY*L*(E*UP**2+G*K**2*U**2)/X70+X6*(-E*IYY 

1  *LAP+E*  IXX*LAP-  IXY*  (G*K**2-E*U**2) )  /  (X5*L)  ♦  (  -G*IYY*K**2*U-G 

2  *IXX*K**2*U)/X2-X12*E*IXY/L**3 

ELSTIF (2, 8)  «  X9*IYY*L*(E*LAP**2+G*K**2*LA**2)/X70+X12*(-E*1XY* 

1  LAP- (E*IYY*U**2+G*IXX*K**2)/X2)/(X6*L) -X12*E*IXX/L**3 

ELSTIF{2.9)  ■  (-X3)*iyY*L*(E*KP*LAP+G*K**3*U)/X20+(-E*IYT*K*LAP 
1  ♦E*IYY*KP*U-G*IXY*K**3)/X2+(E*IYY*K*U+E*IXY*KP)/L 


ELSTIF(2,10)  *  X13=»IYYfL**2*(E*LAP**2+G*K**2*U**2)/X420-(E*IXY* 
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1  LAP+  (E*  IYY*U**2+G*  IXX*K**2)  /X2)  /X6-X6*E*  IXX/L**2 

ELSTIF(2.n)  =  X13*IXY*L**2*(E*LAP*f2+G*K**2*LA**2)/X420+(E*lYY* 

1  UP-E*IXX*LAP+IXY*(G*K»*2-E*U*f2))/X10+L*(G*IYY*K**2*U*G*IXX 

2  *K**2*LA)/X10+{E*IYY*LA+E*IXX*LA)/L+X6*E*IXY/L**2 

ELSTIF(2.12)  *  X3*L*(E*IXY*K»UP-G*IYY*K*U**2)/X2a+(-IXY*(G*K* 

1  U+E*K*LA)-G*IXY*K*U)/X2+(-G*IXX*K-Et'IXX*K)/L 

ELSTIF(3.3)  =  IYY*(E*KP**2+G*Kt-*4)*L/X3*E*(IYY*K**2+A)/L-E*IYY*K 
1  *KP 

ELSTIF(3.4)  =  IYY*L**2t'(E*KP*UP+G*K**3*LA)/X20-L*(E*IYY*K*LAP-E* 
1  IYY*KP*LA+G*IXY*K**3)/X12-E*IXY*K/L+E*IXYfKP 

ELSTIF(3.5)  =  IXY*L**2*(E*KP*LAP+G*K**3*LA)/X20+L*(-E*IXY*K*LAP+E 
1  *IXY*KP*LA+G*IYY*K**3+A*E*K)/X12+E*IYY*K/L'E*IYY*KP 

ELSTIF(3.6)  =  L*(G*IYYfK**2*LA-E*IXY*K*KP)/X3-(G*IXYt‘K**2-E*IXY* 

1  K**2)/X2 

ELSTIF(3,7)  =  X3*IXY*L*(E*KP*LAP+G*K**3*U)/X20+(-E*IXY*K*LAP+E* 

1  IXY»KP*U+G* IYY*K**3+A*E*K) /X2+ (E* IYY*KP-Ef IXY*K*LA) /L 

ELSTIF(3.8)  »  (-X3)*IYY*L*(E*JCP*UP+G*K**3*LA)/X20+(E*IYY*K*UP- 
1  E*IYY*KP*LA+G*IXY*K**3)/X2+<E*IYY*K*LA+E*1XY*ICP)/L 

ELSTIFO.Q)  »  IYY*(E*KP**2+G*K**4)*L/X6“E*(IYY*K**2+A)/L 

ELSTIF(3,10)  «  -IYY*L**2*(E*KP*UP+G*K**3*LA)/X30-L*(-E*IYY*K*LAP 
1  +E*IYY*KP*U-G*IXY*K**3)/X12-*-E*IXY*K/L 

ELSTIF(3.11)  *  -IXY*L**2*(E*KP*UP+G*K**3*LA)/X30+L*(E*IXY*K*LAP- 
1  E*IXY*KP*LA-G*IYY*K**3-A*E*K)/X12-E*IYY*K/L 

ELSTIF(3,12)  =  L*(G*IYY*K**2*U-E*IXY*Kt‘KP)/X6+(G*IXY*K**2+E*IXY* 
1  K**2)/X2 

ELSTIF(4,4)  *  IYY*L  *3*(E*LAP**2+G*K**2*U**2)/X105+X4*L*(E*IXY* 

1  LAP-»(E*IYY*LA**2+G*IXX*K>**2)/X2)/X16*E*IXY*U>»X4*E*IXX/L 

ELSTIF(4,B)  ■  IXY*L**3*(E*LAP**2+G*K**2*U**2)/X106+(-I2)*L*(E* 

1  IYY*LAP-E*IXX*LAP+IXY*(G*K**2-E*U**2))/X1B+(E*1XX*U-E*ITY* 

2  LA)/X2-X4*E*IXY/L 


f 

> 
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ELSTIF(4.6)  =  -L**2*(E*IXY*K*LAP-G*IYY*K*LA**2)/X20-L*(IXY*(G*K* 

1  LA+E*K>i-LA)  +G*  IXY*K*U)  /X12-E*  IXX*K 

ELSTIF(4.7)  =  X13*IXY*L**2*(E*LAP**2+G*K**2*LA**2)/X420-(-E*IYY* 

1  LAP+E*IXX*UP-IXY*(G*K*f2-E*LA**2))/X10-L*(-G*IYY*K**2*LA-G* 

2  IXX*K**2*LA)/X10-(-E*IYY*U-E^IXX*U)/L+X6*E*IXY/L**2 

ELSTIF(4.8)  =  (-X13)*IYY*Lf‘*2*(E*UP**2+G*K**2*LA**2)/X420-(-E* 

1  IXY*UP - (E* I YY*LA  »  *2+G* IXX*K**2) /X2) /X6+X6*E*IXX/L**2 

ELSTIF(4.9)  =  IYY*L**2*(E*KP*LAP+G*K**3*U)/X30-L*(-E*IYY*K*UP+E 
1  *IYY*KP*U-G*IXY*K**3)/X12+E*IXY*K/L 

ELSTIF(4,10)  =  -IYY*L**3*{E*LAP**2+G*K**2*LA**2)/X140+L*(-E*IXY* 

1  LAP- (E*IYY*LA**2+G*IXX*K**2)/X2)/X15+X2*E*1XX/L 

ELSTIF(4.11)  =  -IXY*L*-»3*(E*LAP**2+G-^K**2*U**2)/X140-L*(-E*IYY* 

1  LAP+E*IXX*LAP-IXY*(G*K**2-E*U**2))/X30-L**2*(G*lYY*K**2*U-*-G* 

2  IXX*K**2*LA)/X60-(E»IYY*LA+Ef‘IXX*LA)/X2-X2*E*IXY/L 

ELSTIF(4.12)  =  -L**2-(E*IXY*K*LAP-G*IYY*K*U**2)/X30-L*(-IXY*(G*K 
1  *LA+E*K*LA)-G*IXY*Ki‘LA)/Xl2 

ELSTIF(S,5)  =  L**3*(E*(IXX*LAP**2+A*K**2)+G*IXX*K**2*U**2)/X105+ 

1  X4*L*  (  (E*IXX*U**2•^Q*^Y*K**2)/X2-E*lXY*U?)/XlS'E*Uy*U*X4* 

2  E*IYY/L 

ELSTIF(5,6)  =  L**2*(G*1XY*K*U**2-E*IXX*K*LAP)/X20+L*(G*IYY*K*U- 
1  G‘*'IXX*K*LA-E*IXX*K#U)/X12+E*IXY*K 

ELSTIF(5.7)  =  X13*L*#2*(E*(IXX*LAP**2+A*K**2)+G*IXX*K**2*U**2)/ 

1  X420+(E*IXY*LAP-(E*IXX*LA**2+G*IYY*K**2)/X2)/X5-X6*E*IYY/L**2 

ELSTIF(6.8)  =  (-X13)*IXY*L**2*{E*LAP**2+G*K**2*U**2)/X420+(-E* 

1  IYY*LAP+E*IXX*LAP-IXY*(G*K*'»2-E*LA**2))/X10+L*(G*iyY*K**2*LA 

2  +G*:«X*K**2*U)/X10+(E*IYYt'U+E*IXX*U)/L-X6*E*IXy/L**2 

ELSTIF(8.9)  =  IXY*L**2’*(E*KP*LAP+G*K**3*LA)/X30+L*(E*IXY*K*LAP-E* 

1  IXY*KP*U-G*IYY*K**3-A*E*K)/X12-E*ITY*K/L 

ELSTIF(6.10)  -  -rXY*L**3*(E*LAP**2+C*K**2*LA**2)/X140-L*(-E*IYY* 

1  LAP+E*IXX*LAP-IXY*(G*K**2-E*LA**2))/X30-L**2*(-G*IYY*K**2*U- 

2  G*IXX*K*f2*U)/X60-(-E*IYY*U-E*IXX*U)/X2-X2*E*IXY/L 

ELSTIF(6,11)  =  -L**3*(E*(IXX*LAP**2+A*K**2)+G*IXX*K**2fLA**2)/X140 
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1  +L*(E*IXY*LAP-  (E*IXX*U**2+G*IYY*K*i‘2)/X2)/X15+X2*E*IYY/L 

ELST1F(5.12)  »  L**2*(G*IXYfK*U**2-E*IXX*K*LAP)/X30+L*(-G*IYY*K* 

1  U+G*IXX*K*U-*-E*IXX*K*U)/X12 

ELSTIF(6.6)  =  L*(IXX*(G*U**2+E*K**2)+G*IYY*U**2)/X3+G»(IYY+IXX) 
1  /L 

ELSTIF(6.7)  *  X3*L*(G*IXY*K*U**2-E*1XX*K*LAP)/X20+(G*IYY*K*U-G 
1  *IXX*K*U-E*IXX*K*U)/X2+(-GfIXY*K-E*IXY*K)/L 

ELSTIF(6.8)  =  X3*Lf(E*IXY*K*LAP-G*IYY*K*U**2)/X20+(IXY*{G*K*U+ 

1  E*K*U)  +Gt‘  IXY*K*U)  /X2+  (  -G*IXX*K-E*  IXX*K)  /L 

ELSTIF(6,S)  =  L*(G*IYY*K**2*U-E*IXY*K*KP)/X6+(-G*IXY*K**2-E*IXY* 
1  K**2)/X2 

ELSTIF(6.10)  =  L**2*(E>-IXY*K*LAP-G*IYY*K*U**2)/X30-L*(-IXY*(G*K* 
1  LA+E*K‘!‘LA)-G*IXY*K*U)/X12 

ELSTIF(6.11)  =  L*(-G*IYYt‘K*U-*-G*IXX*K*U+E*IXX*K*U)/X12-L**2*(G* 
1  IXY*K*LAf*2-E*IXX*K>i>LAP)/X30 

ELSTIF(6,12)  =  L*(IXX*(G*U**2+E*K=o*2)+G*IYY*LA**2)/X6-G*(IYY+IXX 
1  )/L 

ELSTIF(7,7)  *  X13*L*(E*(IXX*LAP**2+A*K**2)>*’G*1XX*K**2*U**2)/X36 

1  ♦X12*((E*IXX*U**2+G*IYY*K**2)/X2-E-*IXY*UP)/(X5*L)+E*IXX*U 

2  *LAP+G*IXY*K**2*LA+X12*E*iyY/L**3 

ELSTIF(7.8)  *  (•X13)*IXY*L*(E*UP**2+G*K**2*U**2)/X35+(-X2*E*IXY 

1  *LA*LAP-G*IYY*K**2*U+G*IXX*K**2*LA)/X2*X6*(E*IYY*LAP-E*IXX* 

2  LAP+IXY*(G*K**2-E*U**2))/(X5*L)+X12*E*IXY/L**3 

ELSTIF(7.9)  s  X7*IXY*L*(E’»KP*LAP+G*K**3*U)/X20+(E*IXY*K*LAP+E* 

1  IXY*KP*LA+G*IYY*K**3-A*E*K)/X2+(E*IXY*K*LA-E*iyY*KP)/L 

ELSTIF(7,10)  »  (-Xll)*IXY*L**2*(E*LAP**2+G*K**2*U**2)/X210-(-E* 

1  IYY*LAP+X1 1 *E* IXX*LAP - IXY* (G*K**2-E*LA**2) ) /XIO-L* (G* IYY*K**2 

2  *U*G*IXX*K**2*U)/X10-(E*m*U*E'*IXX*U)/L+X6*E*IXY/L**2 

ELSTIF(7.11)  ■  (-X11)*L**2*(E*(IXX*LAP**2+A*K**2)+G*IXX*K**2* 

1  U**2)/X210+(X«*E*IXY*LAP-(E*IXX*U**2+0*IYY*K**2)/X2)/X5*X6*E 

2  *IYY/L**2 
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ELSTIF(7.12)  «  X7*L*(G*IXYfK*U**2-E*IXX*K*LAP)/X20+(G*iyY*K*U+ 

1  G*IXX*K*LA-E*IXX*K*U)/X2+(G*IXY*K+E*1XY*K)/L 

ELSTIF(8.8)  =  X13*IYY*L*(E*LAP**2+0*K**2*LA**2)/X35+X12*(E*IXY* 

1  UP+(E*1YY*U**2+G*IXX*K**2)/X2)/(X5*L)+E*IYY*LA*LAP-G*IXY* 

2  K**2*U+X12*E*IXX/L**3 

ELST1F(8.0)  *  (-X7)*IYY*L*(E*KP*LAP+G*K**3*U)/X20+(-E*IYY»K*UP 
1  -  E*  I  YY*KP*U+G*  IXY*K**3)  /X2+  (  -E*IYY*K*LA-E*  IXY*KP)  /L 

ELSTIF(8.10)  *  Xll*IYY*L**2*(E*LAP**2fG*K**2*U**2)/X210-(-X6*E* 

1  IXY*UP-(E*IYY*LA**2+G*1XX*K**2)/X2)/XB*X6*E*IXX/L**2 

ELST1F(8.11)  *  Xll*IXY*L**2*(E*LAP**2+G*K**2*LA**2)/X210+(-Xll*E* 

1  IYY*LAP+E* IXXfLAP - IXY* (GfK**2-E*LA**2) ) /XIO+L* ( -G*IYY*K**2*LA- 

2  G*IXX*K**2*LA)/X10+(-E*IYY*U-E*IXX*U)/L-X6*E*IXY/L**2 

ELSTIF(8.12)  *  X7*L*(E*IXY*K*LAP-G*IYY*K>»LA**2)/X20+(IXY*(G*K*U 
1  +E*K*LA)  -G>!‘IXY*K*LA)/X2+(Gi'IXX*K+E*IXX*K)/L 

ELSTIF (9.9)  =  lYY* (E*KP**2+G*K*t'4) *L/X3+E* (IYY*K**2+A) /L+E*IYY*K* 
1  KP 

ELSTIF(9.10)  =  -IYY*L’*«*2*(E*KP*UP+G*K**3*U)/X20-L*(E*IYY*K*UP- 
1  E*IYY*KP*LA+G*IXY*K**3)/X12-E*IXY*K/L-E*IXY*KP 

ELSTIF  (9. 11)  =  -IXY*L*>«‘2*(E*KP*UP+G*K**3*LA)/X20+L*(-E*IXY*K*LAP 
1  +E*rxy*KP*LA+G*IYY*K*’»3+A*E*K)/X12+E*IYY*K/L+E*IYY*KP 

ELSTIF(9.12)  =  L*(G*IYY*K**2*LA-E*IXY*K*KP)/X3+(G*IXY*K**2-E*IXY* 
1  K**2)/X2 

ELSTIF (10. 10)  =  IYy*L**3*(E*LAP**2*G*K**2*U**2)/X106+X4*L*(E* 

1  IXY*LAP+(E*IYY*LA**2+G*IXX*K**2)/X2)/X15-E*IXY*LA+X4*E*IXX/L 

ELSTIFdO.ll)  =  IXY*L*t‘3*(E*LAP-t‘*2+G*K**2*U**2)/Xl05+(-X2)*L*(E 

1  *IYY*LAP-E*IXX*LAP+IXY*(G*K*f2-E*U**2))/X16-(E*IXX*U-E*IYY* 

2  U)/X2-X4*E*IXY/L 

ELSTIF(10,12)  «  L**2*(E*IXY*K*LAP-G*IYY*K*U**2)/X20-L*(IXT*(6*K* 
1  U*E*K*U)+G*IXY*K*U)/X12+E*IXX*K 

ELSTIFdl.ll)  »  L**3*(E*(IXX*LAP**2+A*K**2)+C*IXX*K**2*LA**2)/ 

1  X10B4X4*L*( (E*IXX*U**2+G*IYY*K**2)/X2-E*IXY*LAP)/X16+E*IXY* 

2  U+X4*E*IYY/L 


146 


ELSTIF(11,12)  »  •L**2*(G*IXY*K*U**2“E*IXX*K*LAP)/X20+L*(G*IYY*K* 
1  LA-G*IXX*K*LA-E*IXX-fK*U)/X12-E*IXY*K 

ELSTIF(12.12)  »  L*(IXX*<G*U**2+E*K-*'*2)+G*IYY*U**2)/X3+G*(IYY+ 

X  IXX)/L 
C 

C  ASSIGN  REMAINING  ELEMENTS  IN  SYMMETRIC  FASHION 

C 

DO  300  NI  «  2,  12 
DO  300  NJ  *  1,  NI-1 

ELSTIF(NI.NJ)  »  ELSTIF(NJ.NI) 

300  CONTINUE 
C 

C  COMPUTE  GLOBAL  STIFFNESS  MATRIX  (IN  THE  TANGENT  SYSTEM) 

C 

DO  301  NI  *  1.  12 
DO  301  NJ  *  1.  12 

GLSTIF(I+NI.J+NJ)  =  GLSTIF(I+NI.J+NJ)  +  ELSTIF(NI ,NJ) 

301  CONTIinjE 
C 

C  ASSIGN  VALUES  TO  TG  (TANGENT  TO  GLOBAL  TRANSFORMATION  MATRIX) 

C  EVALUATED  FOR  NODE#  1  OF  THE  ELEMENT 

C 

TG(l.l)  =  -DCOS(PP) 

TG(1.2)  *  S1*DSIN(PP) 

TG(1.3)  »  -C1*DSIN(PP) 

TG(2.1)  *  -DSIN(PP) 

TG(2.2)  »  -S1*DC0S(PP) 

TG(2.3)  =  C1*DC0S(PP) 

TG(3.1)  *  XO 
TG(3.2)  =  Cl 
TG(3.3)  =  SI 
C 

C  COMPUTE  TC  (TANGENT  TO  CHORD  TRANSFORMATION  MATRIX)  FOR 

C  NODE#  1  OF  THE  ELEMENT 

C 

NI  *  3 

CALL  MULMATX(GC,  TG.  TC,  NI,  NI.  NI) 

C 

C  FILL  IN  CORRESPONDING  ELEMENTS  OF  TRANS 

C 

DO  302  NI  >  1.  3 
DO  302  NJ  -  1.  3 

TRANS (I+NI,J+NJ)  »  TRANS (I+NI,J+NJ)  ♦  TC(NI,NJ) 
TRANS(I+NI+3,J+NJ+3)  «  TRANS(I*NI+3, J-»NJ+3)  +  TC(NI.NJ) 
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302 


C 

C 

C 


C 

C 

C 


1 

303 


C 

C 

C 


C 

C 

C 

C 


C 

C 

C 

C 


C 

C 

C 


CONTINUE 

INVERT  TC  TO  OBTAIN  TCINV  FOR  NODE#  1  OF  THE  ELEMENT 

IDGT  =  0 
lA  *  3 
NI  »  3 

CALL  LINVIFCTC,  NI.  lA.  TCINV,  IDGT,  WKAREA,  lER) 

FILL  IN  CORRESPONDING  ELEMENTS  OF  TRANSINV 

DO  303  NI  »  1,  3 
DO  303  NJ  *  1,  3 

TRANSINV(I+NI.J+NJ)  =  TRANSINV(I+NI , J+NJ)  *  TCINV(NI,NJ) 
TRANSINV(I+NI+3,J+NJ+3)  =  TRANSINV(I+NI+3, J+NJ+3) 

+  TCINV(NI.NJ) 

CONTINUE 

COMPUTE  PP  FOR  NODE#  2  OF  THE  ELEMENT 
pp  s  pp  +  P*L/LTUBE 

ASSIGN  VALUES  TO  TG  (TANGENT  TO  GLOBAL  TRANSFORMATION  MATRIX) 
EVALUATED  FOR  NODE#  1  OF  THE  ELEMENT 

TG(l.l)  *  -DCOS(PP) 

TG(1,2)  =  S1*DSIN(PP) 

TG(1,3)  =  -C1*DSIN(PP) 

TG(2.1)  =  -DSIN<PP) 

TG(2.2)  =  -S1*DC0S(PP) 

TG(2.3)  =  C1*DC0S(PP) 

TG(3,1)  =  XO 
TG(3.2)  =  Cl 
TG(3,3)  =  SI 

COMPUTE  TC  (TANGENT  TO  CHORD  TRANSFORMATION  MATRIX)  FOR 
NODE#  2  OF  THE  ELEMENT 

NI  »  3 

CALL  MULMATX(GC.  TG,  TC,  Nl.  NI,  NI) 

FILL  IN  CORRESPONDING  ELEMENTS  OF  TRANS 

DO  304  NI  «  1,  3 
DO  304  NJ  >  1,  3 


I 
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TRAllS(I+NI+6,J+NJ+6)  =  TRANS(I+Nl+6, J+NJ+6)  +  TC(NI.NJ) 
TRANS{I+NI+9.J+NJ+9)  *  TRANS(I+NI+9. J+NJ+9)  +  TC(NI.KJ) 

304  CONTINUE 
C 

C  INVERT  TC  TO  OBTAIN  TCINV  FOR  NODE#  2  OF  THE  ELEMENT 

C 

IDGT  =  0 
lA  »  3 
NI  =  3 

CALL  LIJIV1F(TC.  NI.  lA.  TCINV.  IDGT,  WKAREA,  lER) 

C 

C  FILL  IN  CORRESPONDING  ELEMENTS  OF  TRANSINV 

C 

DO  30S  NI  *  1.  3 
DO  306  NJ  *  1.  3 

TRANSINV(I+NI+6. J+NJ+6)  *  TRANSINV(I+NI+6. J+NJ+6) 

1  +  TCINV(NI.NJ) 

TRANSINV(I+HI+d.J+NJ+9)  =  TRANSINV(I+NI+9. J+NJ+9) 

1  +  TCINV(NI.NJ) 

305  CONTINUE 
C 

C  INCREMENT  CONNECTIVES  I  Af®  J  BY  DOF/NODE  *  6 

C 

1*1+6 
J  *  J  +  6 
C 

C  LOOP  BACK 

C 

200  CONTIinJE 
C 

C  OBTAIN  CHDSTF,  GLOBAL  STIFFNESS  MATRIX  IN  THE  CHORD  SYSTEM.  AS: 
C  CHDSTF  =  TRANS  *  GLSTIF  *  TRANSINV 

C 

NI  a  6*(NELEM+1) 

CALL  MULMATX (TRANS.  GLSTIF.  TEMP24.  NI.  NI.  NI) 

NI  a  6*(NELEM+1) 

CALL  MULMATX (TEMP24,  TRANSINV.  CHDSTF.  NI,  NI.  NI) 

C 

C  OBTAIN  REDFLEX.  A  SUBMATRIX  OF  THE  GLOBAL  FLEXIBILITY 
C  MATRIX  BY  INVERTING  A  SUBMA1RIX  OF  CHDSTF  OF  DIMENSIONS 
C  6aNELEM  x  B+NELEM,  I.E..  WHICH  EXCLUDES  ELEMENTS 
C  CORRESPONDING  TO  THE  FIRST  NODE  WHOSE  DISPUCENENTS  ARE 
C  KNOWN  TO  BE  ZERO 
C 


DO  210  NI  a  7.  6*(NELEM+1) 
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DO  210  NJ  *  7.  6*(NELEM+1) 

TEMP18(NI-6.NJ-6)  *  CHDSTF(NI.NJ) 

210  CONTINUE 
IDGT  =  0 
lA  =  6*NELEM 
NI  =  6+NELEM 

CALL  LINV1F(TEMP18.  NI.  lA.  REDFLEX,  IDGT.  WKAREA,  lER) 

DO  900  NI  =  7.  6*(NELEM+1) 

DO  900  NJ  =  7.  6*(NELEM+1) 

TEMP18(NI-6.NJ-6)  =  CHDSTF(NI.NJ) 

900  CONTINUE 

NI  =  6*NELEM 

CALL  MULMATX(REDFLEX.TEMP18. CHECK. NI.NI.NI) 

C 

C  OBTAIN  DISPLACEMENTS  DELCHD  =  REDFLEX  *  FKNOWN  IN  CHORD  SYSTEM 
C 

NI  =  6*NELEM 
NJ  =  1 

CALL  MULMATX (REDFLEX.  FKNOWN.  DELCHD.  NI.  NI.  NJ) 

C 

C  COMPUTE  REACTIONS  FROM  THE  ABOVE  DISPLACEMENTS 
C 

DO  903  NI  =  1.  6 

DO  903  NJ  *  7.  6*(NELEM+1) 

KBETAALPHA(NI.NJ-6)  »  CHDSTF(NI.NJ) 

903  CONTINUE 
NI  »  6 

NJ  =  6*NELEM 
NK  =  1 

CALL  MULMATX (KBETAALPHA. DELCHD. TEMP61.NI.NJ.NK) 

C 

C  WRITE  COMPUTED  REACTIONS  IN  CHORD  SYSTEM  TO  FILE  "FEMOUTPUT" 

C 

WRITE(12,*) 

WRITE (12.*)  ’COMPUTED  REACTIONS  IN  CHORD  SYSTEM’ 

WRITE(12.*) 

WRITE(12.*)  (TEMP61(NI.l),  NI  *  1.  6) 

WRITE(12.*) 

C 

C  CONVERT  RADIANS  INTO  DEGREES 
C 

DO  212  NI  >  1.  NELEM 
DO  212  NJ  «  1.  3 

DELCHD ((6*NI-NJ+l).l)  «  DELCHD((6*NI-NJ+1) .1)*X180/PI 
212  CONTINUE 
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C 

C  WRITE  CHECK  TO  FILE  "FEMOUTPUT" 

C 

C  WRITE(12.*)  ’FLEX  CHECK:* 

C  DO  901  MI  =  1.  6*MELEM 

C  WRITE(12.*)  (CHECK(NI.NJ).  NJ  =  1.  6*NELEM) 

C  901  CONTINUE 
C 

C  WRITE  END  DISPLACEMENT  TO  FILE  "FEMOUTPUT” 

C 

WRITE(12.*) 

WRITE(12,*)  ’END  DISPLACEMENT  IN  CHORD' SYSTEM: ’ 
WRITE(12.*) 

WRITE(12.*)  ’D«lx  =  ’.  DELCHD(6*(NELEM-1)+1.1).  ’m’ 
WRITE(12.*)  ’D«ly  *  ’.  DELCHD(6*(NELEM-1)+2.1) .  ’m’ 
WRITE(12.*)  ’D«lz  *  ’.  DELCHD(6*(NELEM-1)+3.1).  ’m’ 
WRITE(12.*)  ’Th«x  *  ’,  DELCHD(6*(NELEM-1)+4.1) .  ’d#g’ 
WRITE(12.*)  ’They*  ’.  DELCHD(6*(NELEM-1)+6.1) ,  ’deg’ 
WRITE(12.^)  ’Thez  *  ’.  DELCHD{6*(NELEM-1)+6,1) .  ’deg’ 

C 

C  WRITE  GLOBAL  STIFFNESS  MATRIX  INTO  FILE  "FEMOUTPUT" 

C 

C  WRITE(12.*) 

C  WRITE(12.*)  ’GLSTIF:’ 

C  WRITE(12.*) 

C  DO  214  NI  »  1.  6*(NELEM+1) 

C  WRITE(12.*)  (INT(GLSTIF(NI.NJ)).  NJ  *  1.  6*(NELEM+1)) 

C  214  CONTINUE 
C 

C  WRITE  FKNOWN,  GLOBAL  FORCE  VECTOR  TO  FILE  "FEMOUTPUT" 

C 

WRITE(12.*) 

WRITE(12.*)  ’FKNOWN  VECTOR: ’ 

WRITE(12.*) 

WRITE(12.*)  (FKNOWN(NI.l).  NI  =  1.  6*NELEM) 

C 

STOP 

END 

C 

C  MULMATX  SUBROUTINE  BEGINS 

C 

SUBROUTINE  MULMATX (MATl .MAT2. MATS. LI .L2.L3) 

INTEGER  LI,  L2,  L3,  KNTl,  KNT2.  KNT3 
REAL*8  MAT1(L1,L2),  MAT2(L2,L3}.  MAT3(L1.L3) 

DO  10  KNTl  ■  1.  LI 
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DO  10  K1IT2  =  1.  L3 

MAT3(KNT1.KNT2)  =  O.ODO 
DO  11  KNT3  =  1,  L2 

MAT3(KNT1.KNT2)  *  MAT3(KNT1 ,KNT2)+MAT1(KNT1 ,KNT3)* 
1  MAT2(KNT3.KNT2) 

11  CONTINUE 
10  CONTINUE 
RETURN 
END 
C 

C  END  OF  PROGRAM 

C _ 
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